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PREFACE 


This  progress  report  appears  in  two  parts.  Part  A  is  a  summary  of  work  done  in  support  * 
of  this  program  at  the  Naval  Undersea  Center.  Part  B  contains  final  reports  submitted  by 
contracts  in  support  of  this  program. 

The  major  milestone  achieved  during  this  period  is  the  demonstration  of  the  TV  band¬ 
width  reduction  system  CCD  implementation.  Still  pictures  and  video  tapes  are  produced 
at  various  bit  rates.  System  specifications  are  documented  and  contracts  and  negotiations 
are  in  progress  to  have  a  flight  system  produced.  This  system  is  described  in  Appendix  A. 

Also  accomplished  are  simulation  studies  on  more  advanced  bandwidth  reduction  sys¬ 
tems  that  take  advantage  of  both  spatial  and  temporal  correlation  in  TV  images.  Results 
of  this  work  are  described  in  Appendices  B  and  C. 

The  algorithm  for  the  prime  cosine  transform  is  documented  in  Appendix  D,  and  a  sur¬ 
face  wave  prime  Fourier  transform  device  is  implemented  for  evaluation  of  the  algorithm 
in  Appendix  E. 

Hardware  developed  under  this  program  is  shown  to  have  applications  in  many  areas  of 
signal  processing  and  techniques  for  modularly  expanding  the  basic  building  blocks  of  this 
hardware  into  larger  systems  are  investigated  in  Appendix  F.  Several  signal  processing 
architectures  are  described  in  Appendix  G. 

Simulation  studies  by  USC  and  others  document  that  the  Cosine  transform  is  very  close 
to  the  Karhunen-Loeve  transform  when  applied  to  image  data.  These  results  are  obtained 
from  experimental  data  but  as  yet  no  analytical  explanation  is  available.  However,  NUC 
has  undertaken  to  give  a  theoretical  basis  for  the  performance  of  the  cosine  transform  and 
the  results  are  presented  in  Appendix  H. 

Papers  have  been  presented  at  several  conferences  on  the  BBD  and  CCD  bandwidth 
reduction  systems  and  are  included  as  Appendices  I,  and  J. 

Work  is  also  conducted  on  a  frequency  hopping  modem  and  some  preliminary  results 
are  presented  in  Appendices  K  and  L. 

For  ease  of  content  location,  the  enclosed  appendices  are  numbered  consecutively  in 
Arabic  numerals  beginning  with  the  title  page  of  Appendix  A  and  continuing  to  the  end 
of  Appendix  L. 
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HYBRID  COSINF/DPCM  TV  BANDWIDTH  REDUCTION  HARDWARE 


by  Robert  Means  and  E.  II  Wrench  Jr. 


Following  is  a  technical  description  of  the  bandwidth  compression  system  developed  at 
the  Naval  Undersea  Center.  This  paper  is  an  excerpt  from  the  specifications  for  a  contract 
currently  under  negotiation  to  produce  flight  hardware  for  the  Army  AQUILA  RPV. 
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Section  I.  Bandwidth  Reduction  System  Specifications 


1.1  The  item  required  is  a  system  which  can  reduce  the  data  required  to  trans¬ 
mit  television  images  over  a  digital  communication  modem  to  the  rates  specified 
below.  The  system  performance,  in  terms  of  picture  quality  and  data  rate,  must 
approach  that  observed  in  computer  simulations.  The  system  must  be  capable  of 
accepting  EIA  monochrome  television  signals. 

1.2  Considerable  simulation,  analysis  and  breadboard  experimentation  have  been 
done  on  the  problem  of  optimizing  the  bandwidth  compression  that  can  be  obtained 
by  a  real  time  system  within  a  small  power,  weight,  size  and  cost  constraint. 

The  system  will  be  used  in  the  Aquila  mini  remotely  piloted  vehicle.  Figure 
1.1  is  a  block  diagram  of  the  hardware  required.  It  must  accept  composite  video 
in  the  flight  hardware  (transmitter)  and  return  composite  video  at  the  ground 
station  (receiver).  The  data  rates  out  of  the  source  encoder  are  specified 
below.  Contained  at  the  end  of  this  paper  are  schematics  of  the  subsystem  which 
have  been  implemented  and  tested  at  the  Naval  Undersea  Center.  These  are  provided 
for  informational  purposes.  They  should  be  considered  as  design  guidelines  which 
would  aid  the  contractor  in  meeting  system  specifications. 


SOURCE  ENCODER  SOURCE  DECODER 


-  DISPLAY 


Figure  1.1 


1.3  Each  specification  will  be  expanded  upon  in  a  separate  section  where 
appropriate.  Some  of  these  requirements  are  rigid  and  must  be  adhered  to  by 
the  contractor.  Other  requirements  are  such  that  considerable  design  effort, 
invention  and  ingenuity  must  be  expended  by  the  contractor  to  meet  the  perfor¬ 
mance  specifications.  It  is  realized  by  the  government  that  the  required  sys¬ 
tem  imposes  performance  specifications  that  will  require  novel  techniques  to 
achieve  and  that  there  are  tradeoffs  between  power,  weight,  cost  and  performance. 
Specifications  will  be  found  in  each  section.  A  partial  list  follows 

1.4  Timing 

1)  The  number  of  picture  elements  in  each  horizontal  line  must  not  be 
appreciably  less  than  256. 

2)  The  number  of  lines  per  field  must  be  262.5. 

3)  The  frame  rate  must  be  7.5  frames  per  second  with  the  capability  of 
operating  at  a  reduced  frame  rate  with  a  minimum  of  system  redesign.  Some  design 
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effort  should  be  expended  to  insure  that  the  system  would  not  be  completely 
incompatable  with  a  solid  state  sensor  operating  in  the  snapshot  mode. 

4)  The  video  signal  must  be  transformed  in  horizontal  blocks  of  either 
31  or  32  pixels. 

5)  Video  sample  rate  =4.8  MHz. 

6)  System  clock  frequency  =  9.6  MHz. 

7)  There  must  be  a  sync  output  pulse  every  8  fields. 

1.5  Discrete  Cosine  Transform 

1)  The  horizontal  transform  must  be  a  discrete  cosine  transform. 

2)  Dynamic  range  of  output  coefficients  must  be  2®, 

3)  Peak-to-peak  signal  to  peak-to-peak  noise  of  output  coefficients  must 
be  >  2®. 

4)  Transform  coefficients  must  be  accurate  to  within  one  percent  over 
the  whole  dynamic  range  of  the  transform. 

1.6  Differential  Pulse  Code  Modulator 

1)  The  vertical  transform  must  be  a  differential  pulse  code  modulator 
operating  on  cosine  transform  coefficients. 

2)  Quantizer  must  be  within  the  DPCM  feedback  loop. 

3)  Dynamic  range  at  input  must  be  2®. 

4)  Maximum  error  at  output  with  steady  state  input  signal  must  not 
exceed  one  least  significant  bit  when  quantized. at  6  bits  per  coefficient. 

5)  Quantization  levels  must  be  spaced  for  equal  probability  of  occupancy 
at  all  bit  rates. 

6)  DPCM  must  be  capable  of  variable  bit  assignment  per  coefficient  from 
zero  to  six  bits. 

7)  Output  rate  must  be  continuous  serial  capable  of  switching  between 
200,  400,  800  and  1600  kilobits  per  second  on  command  from  the  modem. 

1.7  Ground  Station 

1)  The  ground  station  must  be  implemented  with  digital  hardware. 

2)  Power  and  weight  are  not  critical  factors  in  the  ground  station. 

1.8  Form  Factors 

1)  The  system  must  fit  on  one  (1)  card,  the  size  of  which  is  approxi¬ 
mately  4"  x  8"  x  1/2".  Exact  form  factors  will  be  specified  after  award  of 
contract. 

2)  The  system  must  have  a  weight  of  less  than  one  pound. 

3)  The  system  must  dissipate  less  then  ten  watts. 

4)  Environemntal  specifications  are  those  specified  for  the  Aquila 
Lockheed  RPV. 


1.9  Administrative 

1)  There  will  be  design  review  meetings  held  every  90  days  at  the  con¬ 
tractors  plant. 

2)  Production  cost  estimates  of  the  system  are  requested  with  a  produc¬ 
tion  run  of  1000  a  year  for  5  years  to  be  used  as  a  basis. 

3)  A  bid  is  requested  on  a  package  consisting  of  two  airborne  subsystems 
and  one  ground  station. 

4)  A  bid  is  requested  on  a  package  consisting  of  3  airborne  subsystems 
and  two  ground  stations. 


5)  Proposal  evaluation  criteria  are  listed  in  order  of  importance. 

a)  weight,  size,  power 

b)  performance  and  design 

c)  production  cost 

d)  contractor  experience  and  personnel 

e)  risks  inherent  in  design 

f)  bid  prices 

6)  Delivery  date  is  9-10  months  after  award  of  contract. 
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Section  II.  Timing  and  Synchronization 


2.1  The  TV  bandwidth  reduction  system  must  operate  synchronously  with  the  TV 
camera  and  the  Harris  modem.  The  clocks  required  to  insure  compatibility  have 
been  determined  and  will  be  supplied  by  Harris.  These  clocks  will  be  available 
to  operate  both  in  the  plane  and  on  the  ground  at  synchronize  time  and  phase 
relationships . 

2.2  The  TV  sync  in  the  plane  will  be  supplied  by  the  camera.  On  the  ground 
the  TV  sync  must  be  reconstructed  by  the  bandwidth  reduction  system  from  a 
timing  pulse  transmitted  from  the  plane  and  the  system  clocks.  The  detailed 
operation  of  the  system  timing  is  described  below. 

2.3  In  the  plane,  the  TV  sync  is  generated  by  a  TV  sync  generator  chip  that  is 
driven  by  a  2.0  MHz  clock.  The  2.0  MHz  clock  is  supplied  by  Harris  and  is 
derived  from  the  system  clock.  (2  MHz  =  4.8  MHz  r  12  x  5.)  Each  TV  line  is 
composed  of  130  of  the  2  MHz  clocks  for  a  line  time  of  65.0  ys.  The  V  and  H 
signals  generated  by  the  sync  chip  are  supplied  to  the  bandwidth  reduction 
system. 

2.4  The  basic  sample  clock  rate  for  the  bandwidth  reduction  system  is  4.8  MHz 
and  is  derived  from  a  9.6  MHz  clock  supplied  by  Harris.  The  TV  line  time  cor¬ 
responds  to  312  clocks  at  4.8  MHz.  Of  these,  256  correspond  to  the  active  TV 
line,  and  56  to  retrace  time.  The  even  field  is  processed  as  though  it  were 
identical  to  the  odd  field,  so  the  vertical  resolution  is  262.5  lines. 

2.5  The  phase  relation  of  the  horizontal  sync  pulse  to  the  4.8  MHz  clock  is 
not  known  but  will  remain  constant  since  the  horizontal  sync  is  derived  from 
the  2  MHz  clock  which  in  turn  is  derived  from  the  4.8  MHz  clock.  Therefore 
each  TV  line  will  start  at  the  same  time  relative  to  the  4.8  MHz  sample  clock. 
This  permits  the  4.8  MHz  clock  to  be  used  to  partition  the  active  TV  lines 
into  256  pixel  blocks.  The  pixels  are  numbered  with  0  on  the  left  and  255  on 
the  right. 

2.6  The  effective  frame  rate  of  the  TV  is  reduced  by  a  factor  of  eight  by  only 
coding  1/8  of  a  TV  line  during  each  horizontal  time  interval  (65  ys) .  Eight 

TV  fields  must  be  swept  out  by  the  camera  before  one  entire  field  has  been 
coded.  During  the  first  field,  pixels  0-31  of  all  lines  will  be  coded  and 
transmitted.  During  the  second  field  pixels  32-63  will  be  coded,  etc. 

2.7  The  blocks  of  32  pixels  form  vertical  bars  in  the  TV  field  are  referred  to 
as  stripes.  These  stripes  are  numbered  0-7  with  stripe  zero  containing  pixels 
0-31.  A.  stripe  counter  is  incremented  by  the  vertical  sync  pulse  obtained 
from  the  camera. 
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2.8  A  pulse  that  exactly  describes  the  beginning  of  stripe  zero  (frame  sync) 
is  used  to  obtain  video  synchronization  between  the  plane  and  the  ground.  The 
frame  sync  pulse  (7.5  Hz  rate)  is  sent  to  the  status  modem  in  the  plane  and  then 
transmitted  to  the  ground.  At  the  ground  the  pulse  is  obtained  from  the  status 
modem  with  its  time  relative  to  the  data  preserved.  The  frame  sync  is  then 
used  to  synchronize  the  pixel  counter,  line  counter,  and  stripe  counter  on  the 
ground  with  those  in  the  plane.  It  is  also  used  to  lock  a  TV  sync  chip  (required 
to  produce  composite  video)  on  the  ground  with  the  one  in  the  camera.  Therefore, 
complete  synchronization  between  air  and  ground  video  is  maintained  by  the  frame 
sync  signal. 

2.9  Another  area  where  synchronization  is  required  is  the  interface  between 

the  TV  bandwidth  compression  system  and  the  Harris  modem.  The  interface  requires 
that  serial  data  be  passed  from  the  TV  system  to  Harris  at  four  fixed  rates:  200, 
400,  800,  and  1600  K  bits/sec.  This  corresponds  to,  respectively,  13,  26,  52,  and 
104  bits/TV  line.  The  output  bits  are  to  be  synchronous  with  a  data  clock,  supplied 
by  Harris  at  the  bit  rate.  The  phase  of  this  clock  relative  to  the  4.8  MHz 
clock  is  not  known  but  will  remain  fixed  at  each  of  the  four  frequencies. 

2.10  The  interface  at  the  ground  between  the  modem  and  the  TV  system  is 
identical  to  the  one  just  described.  The  data  from  the  modem  is  presented 
serially  and  synchronously  with  the  data  clock  from  the  modem.  The  data  rate 
at  the  interface  is  determined  by  the  state  of  2  bits  (x&y)  supplied  to  the 
system  by  the  command  and  control  modem. 
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Section  III.  Discrete  Cosine  Transform 


3.1  There  are  two  discrete  cosine  transforms  which  can  be  defined  for  a  finite 
block  size  data  set.  They  have  been  called  the  odd  discrete  cosine  transform 
(ODCT)  and  the  even  discrete  cosine  transform  (EDCT) .  For  a  complete  discus¬ 
sion  of  the  differences  see  Appendix  A. 

3.2  There  exists  many  possible  algorithms  to  implement  the  discrete  cosine 
transform.  It  is  not  the  purpose  of  this  procurement  to  specify  a  given  al¬ 
gorithm.  However,  the  hardware  implementations  of  each  approach  considered 
by  the  Naval  Undersea  Center  will  be  discussed. 


Chirp  Z  Algorithm 

3.3  Consider  the  odd  cosine  transform.  The  chirp  Z  algorithm  computes  the 
coefficients  of  a  32  point  data  block  by  implementing  the  equation 
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A  block  diagram  of  the  hardware  is  shown  in  figure  3.1.  This  algorithm  was 
implemented  by  the  Naval  Undersea  Center  with  charge  coupled  device  transversal 
filters  developed  by  Texas  Instruments.  It  and  the  rest  of  the  bandwidth  com¬ 
pression  system  has  been  demonstrated  with  real  time  video.  The  complete 
schematics  for  the  system  are  attached  at  the  end  of  this  paper. 
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Figure  3.1.  Block  Diagram  of  Chip  Z  Algorithm  Hardware. 
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3.4  The  hardware  consists  of  three  major  blocks:  The  premultipliers,  the 
transversal  filters,  and  the  postmultipliers.  All  three  blocks  performed 
marginally  at  the  4.8  MHz  video  sample  rate.  The  system  performance  was  close 
to,  but  did  not  equal,  that  of  the  computer  simulations. 

3.5  The  breadboard  system  built  by  NUC  used  an  LSI  chip  developed  by  Texas 
Instruments.  All  the  premultipliers  and  differential  current  integraters  were 
placed  on  the  device.  The  transversal  filters  themselves  performed  adequately 
at  the  sample  rate  of  interest  (4.8  MHz),  but  the  peripheral  devices  -  multi¬ 
pliers  and  differential  current  integrators  -  did  not.  In  the  breadboarded 
system  these  were  implemented  by  separate  modules  or  discrete  components.  It 
was  the  off  chip  modules  which  determined  the  size  of  the  laboratory  bread- 
boarded  system.  In  order  to  meet  the  contract  specifications  it  is  necessary 
to  put  as  many  of  these  peripheral  devices  on  the  chip  as  possible;  i.e.  a 
new  CCD  chip  is  required. 


3.6  Since  the  processor  has  a  complete  line  to  process  32  samples,  it  is 
possible  to  capture  the  video  in  an  analog  delay  line  or  memory  at  high  sample 
rate  (4.8  MHz)  and  read  it  out  into  the  processor  at  a  slow  rate  (1.2  MHz). 

This  option  would  relax  the  speed  requirements  on  the  multipliers,  differential 
current  amplifiers,  and  CCD's.  It  will  also  relax  the  requirements  on  the 
DPCM.  It  could,  however,  add  a  separate  piece  of  hardware  to  the  system. 


3.7  The  premultiplier  should  have  an  accuracy  of  at  least  6  bits.  It  is  pre¬ 
sumed  that  the  video  is  limited  to  6  bit  accuracy,  probably  in  the  display. 

The  transform  coefficients  should  have  a  dynamic  range  of  at  least  eight  bits. 


Prime  transform  algorithm 

3.8  The  prime  transform  algorithm  can  be  implemented  by  the  hardware  shown 
in  figure  2.  The  prime  transform  requires  the  use  of  a  permuting  analog  memory, 
a  transversal  filter  and  a  second  permuting  memory  at  the  output.  In  addition, 
it  requires  a  separate  computation  of  the  D.C.  coefficient.  Also  the  first 
video  sample  must  be  treated  different  than  all  the  rest.  An  analog  memory 
chip  has  been  developed  by  the  Reticon  corporation  under  a  contract  by  the 
Naval  Undersea  Center.  This  system  has  not  been  breadboarded  by  the  Naval 
Undersea  Center  and  the  memory  chip  has  not  been  sufficiently  tested.  However 
the  performance  specifications  claimed  by  Reticon  are  attached  in  the  appendices. 
The  prime  algorithm  for  a  Courier  transform,  is 
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and  R  is  a  primitive  root  of  N. 
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Figure  3.2.  Prime  Transform  Algorithm. 
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3.9  The  calculation  of  the  discrete  cosine  transform  could  be  done  by  straight 
forward  digital  techniques.  Two  methods  of  performing  this  digital  cosine 
transform  have  been  investigated.  The  first  method  is  a  direct  implementation 
of  the  algorithm  while  the  second  involves  a  trigonometric  substitution  to 
reduce  the  complexity  of  the  calculation. 


3.10.  The  odd  discrete  cosine  transform  of  the  form 


has  been  directly  iaplemented  for  N-32.  The  design  uses  a  read-only-memory  to 
store  the  basis  vectors  and  performs  the  multiplication  using  LSI  multipliers. 
(See  Figure  3.S.)  Due  to  the  limitation  of  a  clock  frequency  of  9.6  MHz  and  a 
maximum  of  65  usee  in  which  to  perform  the  calculation,  only  the  first  19  coef' 
ficients  can  be  calculated.  In  order  to  obtain  all  32  coefficients,  it  is 
necessary  to  have  at  least  2  parallel  processors  operating  together. 

3.11  The  transform  can  be  Implemented  without  having  to  do  any  multiplies  by 
the  substitution  of  variable  Fn«cos*1|n.  The  transform  then  becomes 

®k  "  CO*  £0  4  2  [cos  (jot  ♦  Fn)  ♦  cos  (fjjry-  -  Fn)]  k  -  0,1... M-l 
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By  using  this  method,  the  original  transform  requires  only  additions.  (See 
Figure  3.4.)  The  values  of  the  cosine,  arccosine,  and  2trnk/2N-l  are  contained 
in  read-only-memories.  This  arccosine  method  has  been  implemented  for  N=32 
and  requires  about  the  same  number  of  parts  as  the  direct  method.  Since  it  is 
limited  by  the  same  constraints  of  clock  frequency  and  calculation  time,  it 
also  can  calculate  only  the  first  19  of  32  coefficients.  However,  a  unit  has 
been  built  using  two  processors  in  parallel  that  can  calculate  all  32  coeffi¬ 
cients  by  having  one  processor  calculate  G0  thru  G15  while  the  second  processor 
calculates  G^  thru  G31  at  the  same  time.  An  overflow  limiting  feature  has 
been  included  in  this  design.  This  feature  makes  it  possible  to  experiment 
with  the  selection  of  bits  to  be  used  for  the  output  without  suffering  large 
errors  from  an  overflow  condition.  This  dual  processor  unit  requires  about  60% 
more  parts  than  the  single  unit.  Both  the  direct  implementation  and  the 
arccosine  method  units  use  double  buffering  at  both  the  input  and  output  so 
that  the  processor  can  be  operating  on  data  while  other  data  is  read  in  and 
results  read  out.  This  increases  the  time  available  for  calculation.  Sche¬ 
matics  for  these  systems  are  attached  at  the  end  of  this  paper. 
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Section  IV.  Prewhitening  Filter 


4.1  There  are  two  reasons  to  use  a  prewhitening  filter  as  the  first  operation 
on  the  video  signal. 

4.2  One:  Since  the  cosine  transform  compacts  the  sampled  video  into  the  low 
order  transform  coefficients,  and  the  subsystem  which  calculates  the  cosine 
transform  has  a  limited  dynamic  range,  the  high  frequency  coefficients  will 
have  a  poor  signal  to  noise  ratio.  One  method  of  increasing  the  signal  to 
noise  ratio  for  the  high  frequency  coefficients  is  to  filter  the  video  to  empha¬ 
sise  the  high  frequency  information.  The  phase  requirements  and  invertability 
of  this  filter  and  their  implications  will  be  discussed  later. 

4.3  Two:  The  DPCM  requires  a  non-linear  quantizer.  The  choice  of  the  quanti¬ 
zation  levels  is  extremely  important  to  the  resulting  picture  quality.  The 
levels  are  chosen  by  examining  the  probability  distribution  for  the  difference 
between  transform  coefficients.  Since  the  variance  of  the  transform  coeffici¬ 
ents  decreases  as  a  function  of  n (coefficient  number),  each  coefficient  requires 
a  separate  quantizer.  This  requirement  can  be  eliminated  by  normalizing  all 
the  coefficients  to  have  the  same  variance.  The  normalization  of  the  coeffici¬ 
ents  can  be  done  by  multiplying  each  coefficient  by  the  appropriate  scale  factor 
as  it  enters  the  DPCM  subsystem.  However,  the  high  frequency  coefficients 
would  still  have  a  poor  signal  to  noise  ratio.  Instead,  since  this  operation 

on  the  transform  coefficients  is  equivalent  to  a  frequency  domain  multiply,  one 
can  inclement  the  scaling  necessary  for  the  DPCM  by  a  prewhitening  filter  direc¬ 
tly  on  the  video. 

Design  considerations  of  the  prewhitening  filter 

4.4  Let  us  assume  that  the  variance  of  the  cosine  transform  coefficients 
decreases  proportional  to  2/(n+2).  (This  is  a  fair  approximation  for  most  TV 
scenes.)  The  "ideal"  prewhitening  filter  would  multiply  the  coefficients  by 
(n«-2)/2.  This  filter  is  very  difficult  to  implement  in  the  time  domain.  How¬ 
ever,  one  can  implement  a  simple  filter  which  approximates  the  "ideal"  function. 
The  difference  between  the  approximation  and  the  "ideal"  equalization  can  be 
removed  by  the  use  of  a  multiplier  before  the  DPCM.  This  multiplier  then  oper¬ 
ates  on  coefficients  with  good  signal  to  noise  ratio  and  only  has  to  multiply 
coefficients  by  a  relatively  small  factor.  A  simple  example  of  a  prefilter 
which  can  be  implemented  by  a  delay  line  and  a  differential  amplifier  is  shown 
in  figure  4.1 

-A7  **- 
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Figure  4.1.  Simple  Example  of  a  Prefilter. 


The  impulse  response  of  this  filter  is  also  shown.  In  the  frequency  domain 
this  filter  multiplies  the  Fourier  transform  coefficients  by  the  transfer  func¬ 
tion 


H (f )  =  1  -a+  2i  sin  irfAt 


4.5  There  are  severe  limitations  on  the  filter  phase  distortion  allowed.  The 
filter  in  the  example  has  complex  phase  shift.  It  is  not  possible  to  compute 
a  cosine  transform  of  the  video  and  use  this  prefilter.  If  ot=1.0,  then  it  is 
possible  to  compensate  for  the  purely  imaginary  phase  shift,  but  it  introduces 
problems  with  dc  signals.  The  filter  should  multiply  the  variance  of  the  cosine 
transform  coefficients  by  the  transfer  function  of  the  filter.  Since  the  cosine 
transform  of  the  video  is  the  Fourier  transform  of  the  extended  data  set  (Appen¬ 
dix  A)  one  can  define: 

g(t)  =  sampled  video  data  set 

ge(t)  =  extended  data  set  =  g(t)  +  g(-t) 

h(t)  =  impulse  response  of  the  prefilter 

F(g(t))  =  Fourier  transform  of  g(t)  =  G(f) 

C(g(t))  =  cosine  transform  of  g(t)  =  F(g(t)  +  g(-t))  =  G(f)  +  G*(f) 


Thus 


C(h(t)  *  g(t))  =  F(h(t)  *  g(t)  +  h(-t)  *  g(-t))  (1) 

=  H(f)  G  (f)  +  H*  (f)  G*  (f) 

=  Re  H(f )  G  (f)  =  Re  F(h  *  g) 


If  H(f)  is  purely  real,  then  the  cosine  transform  of  h  *  g  can  be  expressed  as 


C(h  *  g)  =  H(f)  [G  +  G*] 


(2) 


Thus  we  have  multiplied  the  cosine  transform  coefficients  by  the  transfer 
function  of  the  prefilter.  It  appears  that  a  prefilter  with  a  real  transfer 
function  is  the  simplest  one  to  achieve  the  proper  whitening.  Such  a  filter 
is  slightly  more  complicated  to  build.  The  impulse  response  of  such  a  filter 
is  given  in  figure  4.2 
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Figure  4.2.  Filter  Impulse  Response. 


In  the  frequency  domain  this  filter  has  the  transfer  function 


H(f)  =  l-2a  cos  2irfAt 

4.6  It  has  been  assumed  that  the  inverse  filter  will  be  implemented  digitally 
in  the  ground  station.  The  inverse  of  the  filter  shown  in  figure  4  cannot 
be  implemented  easily.  However,  this  is  at  the  ground  station.  The  impulse 
response  is  shown  in  figure  4.3 

1.0  , 


Figure  4.3 


where 

3  ■  ^  (1-  /M?)  for  a  <  0.5 


This  inverse  can  be  approximately  implemented  by  the  hardware  shown  in  figure  4.4 

£ 


4.7  The  accuracy  of  this  filter  will  depend  on  the  njmber  of  stages  in  the 
nonrecursive  part  and  the  value  of  6.  For  example:  for  a  pre-emphasis  filter 
with  a*0.4,  then  B*C.5  and  with  eight  stages  in  the  non-recursive  filter,  the 
accuracy  of  the  output  will  be  greater  than  one  part  in  2®.  Output  is  valid 
after  n+1  stages. 

4.8  It  appears  to  be  preferable  to  implement  the  prefilter  which  has  a  real 
transfer  function  for  several  reasons:  the  cosine  transform  does  not  have  to 
be  modified;  noise  does  not  accumulate  in  the  inverse  filter;  and,  most  impor¬ 
tant,  a  dc  bias  on  the  video  signal  reaches  the  cosine  transform  as  dc. 
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4.9  The  discussion  of  the  prefilter  in  this  section  has  gone  into  detail  on 
a  number  of  alternative  designs.  Attention  was  paid  to  the  impact  on  the 
algorithms.  It  was  not  meant  to  give  the  impression  that  the  prefilter  design 
is  extremely  critical  to  picture  quality.  Several  different  prefilters  have 
been  implemented  by  NUC  and  the  differences  in  picture  quality  were  observed 
to  be  small.  It  was  important,  however,  to  implement  some  type  of  pre-emphasis 
and  de-emphasis.  The  pre-emphasis  should  process  the  whole  line  of  video  so 
that  the  de-emphasis  filter  operates  across  the  boundaries  of  the  horizontal 
interlace  used  in  the  system. 


Section  V. 


Differential  Pulse  Code  Modulator 


5.1  The  differential  pulse  code  modulator  (DPCM)  is  the  system  element  that 
produces  the  actual  bandwidth  reduction.  The  processes  of  prewhitening  and 
cosine  transfering  are  exact  and  invertable  and  therefore  no  bandwidth  reduc¬ 
tion  occurs.  The  DPCM  however  is  not  exactly  invertable  because  approximations 
are  made  in  the  quantization  of  the  differences  between  coefficients.  The  use 
of  the  approximate  difference  results  in  a  reduction  in  the  amount  of  data 
required  to  transmit  the  coefficients.  The  transmitted  data  rate  is  varied 

by  changing  the  number  of  quantization  levels  used  for  each  coefficient. 

5.2  The  basic  structure  of  the  DPCM  is  as  shown  below. 


The  operation  is  as  follows.  The  DCT  coefficients  (Gyj' s)  enter  in  time  sequence, 
and  are  differenced  with  a  fraction  (a)  of  the  corresponding  coefficient  from 
the  previous  line  (G^'s).  The  difference  Dn  is  then  approximated  by  one  of  2k 
possible  values  by  the  quantizer.  The  value  of  k  for  each  coefficient  is  pro¬ 
grammed  in  a  ROM.  The  k  bits  from  the  quantizer  that  define  the  quantization 
level  are  then  transmitted  serially  to  the  modem.  The  quantized  difference, 

‘bn,  is  then  added  to  a  •  G_"  to  produce  an  approximation,  Gn ,  of  the  input 
coefficient  (U.  This  valueis  stored  in  the  delay  memory  to  be  used  as  in 
the  next  line. 

5.3  It  is  important  to  note  that  by  placing  the  quantizer  inside  the  feedback 
loop  the  quantization  error  is  kept  from  accumulating,  and  the  error  converges 
to  the  smallest  quantization  level  in  steady  state. 

5.4  The  most  critical  portion  of  the  DPCM  is  the  design  of  the  quantizer.  For 
the  coding  of  the  coefficient  differences  to  be  efficient,  all  quantization 
levels  must  be  equally  likely.  It  has  been  shown  by  Habibi  (Appendix  .1  references) 
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and  others  that  for  a  wide  variety  of  pictures,  the  amplitudes  of  the  coeffi¬ 
cient  differences  are  approximately  exponentially  distributed.  The  probability 
distribution  for  the  differences  is  therefore  assumed  to  be 


f(Dn)  =  an  *  e-CSJDnl 

where  is  the  variance  of  the  nth  coefficient.  Since  the  coefficients  have 
been  prewhitened,  all  c*n's  are  approximately  equal.  The  distribution  is  sym¬ 
metric  about  zero,  so  only  the  positive  quantization  levels  need  be  calculated 
The  probabilities  of  a  difference  falling  into  any  quantization  level  must  be 
equal,  i.e. 

/-bi.l  j 

J tj  f(x)  dx  ’  pr 

k 

where  2  is  the  number  of  quantization  levels  and  bj  is  the  boundary  between 
the  i-1  and  ith  quantization  level.  For  the  positive  quantization  levels 


f(x)  dx 


i 


1  2^"^ 

or  b  -  i  In  [4^— )  (A) 

i  2  -i 

There  are  also  boundaries  at  zero  and  at  negative  b^. 

5.5  Once  the  quantizer  boundaries  are  determined,  the  quantization  values,  Q^,  are 
chosen  to  minimize  the  mean  square  quantization  error,  i.e., 


=  Minimum  where  x^  is  an  input  to  the 
quantizer  that  falls  between  bj_j 
and  bi- 


The  value  of  is  determined  by  setting 
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x.  =  Q.  where  x.  is  the  mean  of  x.  for  contained  ii 
1  1  between  and  in  quantitation  level  i. 


p 

J  b.  xf(x)dx 
Q.  =  — - 

1  f 1 

J  f(x)  dx 
Di-1 


5.6  The  DPCM  should  have  an  input  dynamic  range  of  +2  .  This  implies  that 
the  difference  between  coefficients  has  a  dynamic  range  of  ±2®.  The  input 
gain  of  the  DPCM  should  be  set  such  that  the  probability  of  the  magnitude  of 
the  coefficient  difference  exceeding  2®  is  1/2P  where  p=kffiax  (the  largest 
number  of  bits  that  can  be  used  in  the  quantizer)  i.e. 


-  2® 

L  fWdX*),8  f(X)dX  '  7 


,8  f<X)dX  ' 


-a*2  1 


a  »  2'8  •  In  2?*1 


For  the  DPCM  required  in  this  specification,  p=6.  Thus  a=0.0162.  The  quanti¬ 
zation  boundaries  and  mean  values  can  be  determined  from  equation  A  and  B.  The 
values  are  given  in  Table  1.  It  should  be  noted  that  to  determine  between 
which  quantization  level  a  difference  falls,  for  k  bits  of  quantization  (k<6) , 
it  is  necessary  only  to  truncate  the  six  bit  quantization  number  to  k  bits. 

This  will  become  clear  by  examining  Table  1. 


QUANTIZER 

POSITIVE  VALUES 

FOR  64  LEVELS 

QUANTIZER 

QUANTIZER 

TRANSMITTED 

BOUND 

VALUE 

SITS 

0 

1 

000000 

81 

01 01 11 

2 

85 

3 

000001 

89 

011000 

4 

94 

5 

000010 

98 

011001 

6 

103 

7 

000011 

108 

oiioio 

8 

‘  14 

9 

000100 

121 

011011 

10 

128 

11 

000101 

137 

011100 

13 

146 

14 

000110 

158 

011101 

15 

171 

16 

000111 

190 

011110 

18 

213 

19 

001000 

275 

01 1 1 1 1 

20 

21 

001001 

23 

24 

001010 

28 

27 

ooi  on 

29 

30 

001 100 

32 

33 

001101 

35 

37 

© 

«4 

14 

«4 

o 

o 

39 

41 

001111 

43 

45 

010000 

47 

49 

010001 

51 

53 

010010 

55 

57 

010011 

60 

63 

010100 

66 

69 

010101 

72 

75 

010110 

79 

<* 

TABLE  1.  DPCM  QUANTIZER  VALUES 
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QURM7IZEP  POSITIVE  VALUES  FOP  33  LEVELS 


QUANTIZE# 

QUANTIZE# 

TRANSMIT 

BOUND 

VALUE 

PITS 

0 

3 

ooooo 

4 

o 

6 

00001 

10 

00010 

i 

15 

0001 1 

13 

30 

00100 

33 

36 

00101 

39 

•33 

001 10 

35 

39 

00111 

43 

47 

01 000 

51 

55 

oiooi’ 

60 

66 

01010 

73 

78 

oioii 

35 

94 

01100 

103 

1 15 

oiioi 

133 

147 

OHIO 

171 

333 

01111 

TABLE  1.  DPCM  QUANTIZER  VALUES  (Cont'd) 


QUANT  I  ZER 

POSITIVE  VALUES 

FOP  16  LEVELS 

QUANTIZER 

QUANTIZER 

TRANSMITTED 

BOUND 

VALUE 

BITS 

0 

4 

0000 

8 

13 

0001 

18 

23 

0010 

29 

36 

0011 

43 

51 

0100 

60 

72 

0101 

85 

104 

0110 

128 

190 

Dili 

QUANT I ZER 

POSITIVE  VALUES  FOR  8  LEVELS 

QUANTIZER 

QUANTIZER 

TRANSMITTED 

BOUND 

VALUE 

BITS 

0 

9 

000 

►  » 

00 

30 

001 

43 

62 

010 

85 

147 

011 

TABLE  1. 

DPCM  QUANTIZER  VALUES 

(Cont'd) . 
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QUANTIZER  POSITIVE  VALUED  POP  4  LEVELS 


QUANTIZER 

BOUND 


QUANT  I ZEP 
VALUE 


'PANSMITTED 

BITS- 


QUANTIZEP  POSITIVE  VALUES  POP  £  LEVELS 


QUANTIZER 

bound 


QUANTIZER 

VALUE 


TRANSAIT^ED 

BITS 


TABLE  1.  DPCM  QUANTIZER  VALUES  (Cont'd) 
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5.7  The  inverse  DPCM  at  the  ground  station  is  almost  identical  to  the  DPCM 
only  less  hardware  is  required.  The  inverse  DPCM  is  shown  below. 


5.8  Two  different  designs  for  the  DPCM  have  been  built  at  NUC.  The  first  is 
a  hybrid  analog,  digital  system  shown  schematically  below 


Coc  ( 
N 


-  f- 
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This  design  has  the  advantage  that  the  analog  to  digital  converter  is  nonlinear. 
It  has  8  bits  of  dynamic  range  but  only  4  bits  of  resolution.  At  low  bit  rates 
(i.e.  large  bandwidth  reduction)  4  bits  or  16  quantization  levels  are  all  that 
are  required.  Therefore,  a  very  fast  A/D  and  DPCM  can  be  implemented  using  16 
comparators.  This  DPCM  design  has  two  disadvantages.  1)  The  analog  differ¬ 
ential  amplifier  and  D/A  converter  at  the  input  can  introduce  noise  into  the 
system.  2)  For  high  data  rates  (low  compression)  it  is  desriable  to  have 
more  than  16  quantization  levels.  At  1.6  m  bits/sec,  6  bits  or  64  levels 
are  desirable  for  the  first  few  low  frequency  coefficients. 


5.9  A  second  all  digital  DPCM  has  been  designed  and  is  currently  being  imple¬ 
mented.  It  is  intended  to  be  used  with  a  cosine  transform  that  reduces  the 
processing  speed  from  4.8  Wz  to  1.2  MHz.  A  small  relatively  low  speed  8  bit 
A/D  converter  can  ther  be  used  directly  on  the  output  of  the  DCT.  All  DPCM 
arithmetic  operations  are  then  performed  digitally.  This  approach  uses  some¬ 
what  less  hardware  since  1)  it  is  running  slower  and  2)  only  one  A/D  is  required 
rather  than  both  an  A/D  and  D/A.  The  disadvantage  of  the  design  is  that  the  coeffi¬ 
cients  must  enter  the  DPCM  at  a  low  rate  (1.2  tffz)  or  a  4.8  (Vz  A/D  must  be  used. 
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5.10  The  block  diagram  of  the  all  digital  DPCM  is  shown  below. 


In  both  designs  the  quantization  occurs  in  two  steps.  1)  The  input  data  is 
assigned  to  one  of  2  max  quantization  levels.  2)  The  k  bits  are  then 
truncated  to  k  bits  as  dictated  by  the  ROM  that  stores  ?1ie  number  of  bits  for 
each  coefficient.  These  k  bits  and  the  bits  for  the  nu.  r  of  bits  per  coeffi¬ 
cient  address  a  ROM  that  outputs  the  appropriate  quantization  value.  The  data 
for  the  ROM's  for  a  6  bit  quantizer  are  given  in  Table  1.  The  input  data  rate 
for  the  DPCM  is  determined  by  the  DCT.  It  is  desirable  that  the  DCT  have  an 
output  rate  less  than  the  4.8  MHz  to  reduce  the  speed  at  which  the  DPCM  must 
operate.  The  output  data  rate  for  the  DPCM  is  fixed  by  the  modem  interface 
requirements.  The  DPCM  output  is  serial  at  four  fixed  rates:  200,  400,  800 
and  1600  MHz.  The  required  data  rates  translate  to,  respectively,  13,  26,  52, 
and  104  bits  per  TV  line.  The  appropriate  data  rate  is  obtained  by  choosing  the 
number  of  bits  used  to  quantize  each  coefficient  such  that  the  sum  equals  the 
required  number  of  bits  per  TV  line.  The  optimum  bit  assignments  for  minimizing 
mean  square  error  have  been  determined  by  Habibi  at  U.S.C.  and  are  given  in  the 
table  below.  These  bit  assignment  are  not  necessarily  best,  and  the  DPCM  design 
should  allow  bit  assignment  changes  to  be  made  easily  (i.e.  changing  a  ROM). 

The  DPCM  must  be  able  to  switch  between  the  four  bit  rates  on  command  from  the 
ground.  Two  input  bits,  x  and  y,  are  supplied  by  the  command  and  control  link 
that  specify  the  bit  rate,  as  shown  below. 


x _ £  Bit  Rate  (kilobits/sec) 

0  0  1600 

0  1  800 

1  0  400 

1  1  200 

5.11  The  output  data  from  the  DPCM  is  serial  and  the  bits  are  presented  to 
the  modem  syncronously  with  a  data  clock  at  the  appropriate  rate  (200  KHz, 

400  KHz,  etc.).  The  data  clock  is  supplied  by  the  modem.  The  data  is  generated 
by  the  DPCM  at  a  nonuniform  rate  and  a  small  buffer  is  required  between  the 
DPCM  and  the  modem. 
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5.12  The  interface  from  the  DPCM  to  the  modem  is  currently  implemented  and  is 
described  below.  (Sea  schematics  at  the  end  of  this  report.) 


5.13  As  mentioned  earlier,  when  fewer  than  6  bits  are  used  to  quantize  a 
coefficient,  only  the  high  order  bits  of  the  6  bit  word  are  retained.  The 
interface  is  a  double  buffered  arrangement  where  all  thirty- two,  6  bit  quan¬ 
tization  levels  are  stored.  Only  selected  bits  are  read  from  the  eemory  for 
transmission  to  the  modem.  A  read  only  memory  is  addressed  sequentially  by 
the  data  clock  from  the  modem.  A  ROM  then  reads  the  appropriate  bits  from 
one  buffer,  while  new  data  is  stored  in  the  other  buffer.  At  the  end  of  each 
TV  line  the  rolls  of  the  buffers  are  interchanged.  The  bits  that  should  be 
truncated  because  fewer  than  6  bits  are  used  for  a  coefficient  are  not  addressed 
by  the  R0M~an<T therefore'  not  transmftfedi  The  "order  in  which"  the  RCH  addresses 
the  bits  is  such  that  all  the  bits  that  must  be  sent  at  the  200  KHz  data  rate 
are  transmitted  first.  Next,  all  bits  required  for  the  400  KHz  data  rate  (but 
not  for  the  200  KHz  rate)  are  transmitted.  This_pattern  continues  until  all 
bits  have  been  accounted  for.  Since  the  ROM  address  is  incremented  by  the 

data  clock  from  the  modem  and  reset  to  zero  at  the  beginning  of  each  TV  line, 
the  number  of  bits  sent  is  determined  only  by  the  data  clock  frequency,  and  the 
correct  bits  are  always  transmitted. 

5.14  The  Interface  at  the  ground  station  is  the  reverse  of  the  above  process. 

A  ROM  stores  the  buffer  address  in  which  each  input  bit  is  to  be  stored.  The 
quantized  data  is  then  reconstructed  in  one  buffer  while  the  other  buffer  feeds 
the  inverse  DPCM.  Again,  the  role  of  the  two  buffers  is  changed  at  the  begin- 
nlng  of  each  TV  llne.~ 

5.15  This  represents  only  one  of  many  possible  interfaces  that  could  be  designed 
and  is  probably  not  a  minimum  hardware  implementation.  An  alternative  approach 
is  to  use  a  first-in-first-out  memory  (FIFO)  to  buffer  the  data  between  the 

DPCM  and  the  modem.  This  requires  an  extremely  fast  FIFO  but  could  significantly 
reduce  parts  count. 


■j.16  An  all  digital  DPCM  is  the  more  desirable  implementation  because  it 
allows  for  6  bit  quantization  with  minimum  hardware.  However,  it  does  need 
either  a  fast  A/D  or  pcrfcrably  a  slow  input  rate  from  the  PCT .  The  drawings 
for  a  possible  digital  implementation  that  has  not  been  fully  tested  are  given 
at  the  end  of  this  report.  Schematics  for  the  hydrid  analog/digital  DPCM  are 
available  from  NUC  if  desired. 


Section  VI.  Frame  Store  Memory 


6.1  A  frame  store  memory  (FSM)  is  required  in  the  ground  station  to  reconstruct 
the  picture  from  the  8  stripes  that  are  transmitted  per  frame.  The  frame  store 
memory  as  currently  implemented  uses  a  commercial  32k  x  18  bit  solid  state 
memory  (specs  in  appendix)  and  stores  the  pixels  sequentally  3  pixels  per  18 
bit  word.  All  reconstructed  picture  elements  are  buffered  in  the  FSM  and  read 
out  at  EIA  standard  TV  rates.  Because  of  the  8  to  1  reduction  in  frame  rate, 
each  picture  element  is  displayed  8  times  before  being  updated. 

6.2  The  basic  operation  of  the  frame  store  memory  is  as  follows.  The  TV  sync 

in  the  plane  and  on  the  ground  are  locked  together  as  discussed  in  the  timing  and 
synchronization  section.  During  the  active  part  of  each  TV  line,  32  picture 
elements  (pixels)  are  reconstructed  at  the  ground  and  buffered.  Simultaneously, 
256  pixels  are  read  from  the  main  memory  and  displayed.  During  the  horizontal 
retrace  time,  the  32  new  pixels  are  stored  in  the  main  memory.  Since  the  camera 
and  the  display  are  sync-locked,  the  new  pixels  are  stored  in  the  same  line  as 
was  just  displayed.  The  stripe  (1/8  line)  in  which  the  pixels  are  stored  is 
determined  by  the  stripe  counter  on  the  ground,  which  is  also  synchronous  with 
the  one  in  the  plane. 

6.3  There  is  no  distinction  made  anywhere  in  the  bandwidth  compression  system 
between  the  even  and  odd  fields  of  the  TV  frame.  The  even  and  odd  fields  are 
assumed  to  be  the  same,  and  the  frame  store  memory  stores  only  one  field,  i.e. 

256  pixels  x  262.5  lines.  It  is  possible  by  doubling  the  FSM  size  to  store 
both  fields  and  thus  increase  the  number  of  vertical  lines  to  525.  There  is 
currently,  however,  no  desire  to  do  this. 

6.4  There  are  complications  to  the  operation  of  the  frame  store  memory  due  to 
the  fact  that  3  pixels  are  stored  in  each  main  memory  word.  This  means  that  32 
pixels  require  10  2/3  words.  Since  the  line  is  divided  into  eight  32-pixel- 
stripes,  some  stripes  do  not  start  or  end  on  even  word  boundaries  in  the  main 
memory.  A  problem  arises  when  only  part  of  a  main  memory  word  is  to  be  updated 
(i.e.  6  or  9  bits).  Because  that  part  of  the  word  not  updated  is  destroyed, 
one  must  compensate  by  reading  from  the  main  memory  the  pixels  that  could  be 
destroyed  in  the  update  process  (i.e.  the  first  and  last  words  of  the  block). 

The  new  pixels  are  merged  with  the  data  from  the  memory  and  rewritten,  thus 
saving  the  pixels  that  are  not  updated. 

6.5  The  current  implementation  of  the  FSM  has  only  256  lines  per  frame  and  is 
not  set  up  to  accept  262.5  lines  per  field  video.  The  modification  to  the 
existing  design  to  incorporate  the  262.5  line  video  should  be  very  minor  and 
NUC  can  provide  guidance  in  this  area. 
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COMBINED  SPATIAL  AND  TEMPORAL  CODING 
OF  DIGITAL  IMAGE  SEQUENCES 
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Abstract 

i 

Interframe  coding  of  television  images  encompasses  techniques  which  make  use  of  corre¬ 
lations  between  pixel  amplitudes  in  successive  frames.  Intraframe  coding  techniques  that 
exploit  spatial  correlations  can,  in  principle,  be  extended  to  include  correlations  in  the 
temporal  domain. 

In  this  paper,  successive  frames  of  digital  images  are  coded  using  two-dimensional  spa¬ 
tial  transforms  combined  with  DPCM  in  the  temporal  domain.  Specific  transform  techniques 
investigated  are  the  two-dimensional  cosine  and  Fourier  transforms.  Due  to  DPCM  encoding 
in  the  temporal  domain,  the  hybrid  transform/DPCM  encoders  require  storage  of  only  the 
single  previous  frame  of  data. 

Hardware  implementation  of  the  Fourier  transform  involves  manipulation  of  complex  num¬ 
bers  where  the  cosine  transform  does  not.  However,  the  Fourier  transform  is  attractive 
because  frame-to-frame  motion  compensation  can  be  introduced  directly  in  the  phase  plane 
by  application  of  appropriate  phase  correction  factors. 

Results  are  presented  in  terms  of  coding  efficiency,  storage  requirements,  computational 
complexity,  and  sensitivity  to  channel  noise. 

Introduction 

In  the  design  of  image  coding  systems  for  digital  communications  channels,  the  primary 
objective  is  to  minimize  the  number  of  code  bits  required  to  reconstruct  the  image  at  the 
receiver.  Reduction  in  the  number  of  code  bits  transmitted  results  in  reduced  channel 
bandwidth,  more  rapid  transmission  of  digital  images,  and  lower  transmitter  power  require¬ 
ments  . 

efficient  coding  of  the  digital  images  is  accomplished  by  removal  of  statistical  redun¬ 
dancies  that  exist  within  the  image.  Transform,  predictive,  and  hybrid  transform/predic¬ 
tive  image  coding  techniques  have  been  developed  to  exploit  intraframe  spatial  image  redun¬ 
dancies.  This  paper  describes  efforts  to  extend  these  image  coding  techniques  to  coding  of 
time-sequences  of  digital  images  transmitted  over  a  digital  communications  channel.  The 
emphasis  has  been  directed  towards  definition  of  an  image  coding  system  that  exploits 
temporal  as  well  as  spatial  image  redundancies. 


Intraframe  Image  Coding 


The  primary  techniques  that  have  been  developed  for  intraframe  image  coding  in  the  spa¬ 
tial  domain  are  transform,  linear  predictive,  and  hybrid  transform/linear  predictive  tech¬ 
niques.  Operational  descriptions  of  these  coding  methods  are  given  below. 


Transform  Image  Coding 


The  basic  premise  of  transform  image  coding  is  that  the  transform  domain  representation 
of  an  image  has  an  energy  distribution  that  is  more  compact  and  therefore  easier  to  effi¬ 
ciently  code  than  the  spatial  domain  version.  In  transform  coding  systems,  a  one-  or  two- 
dimensional  linear  transform  of  an  image  line  segment  or  block  is  performed  at  the  coder. 
The  transform  coefficient  statistics  are  computed  prior  to  quantization  and  coding  for 
transmission.  After  decoding  at  the  receiver,  an  inverse  transform  it.  taken  to  obtain  a 
reconstructed  image.  Transforms  that  have  proven  useful  for  this  application  include 
Fourier,  cosine,  Hadamard,  Slant,  and  Karhunen-Loeve. 11,2,3]  it  should  be  noted  that  the 
two  transforms  can  be  different  for  the  two  spatial  directions. 


Two-dimensional  transforms  have  the  inherent  disadvantage  of  requiring  an  intermediate 
memory  to  store  the  transform  coefficients  computed  in  one  direction  while  the  transform 
is  being  computed  in  the  other  direction.  Transform  coding  techniques  have  been  explored 
extensively  both  theoretically  and  by  simulation.  It  has  been  shown  that  significant  bit 
rate  reductions  can  be  achieved  in  many  applications  with  minimal  image  degradation. 


•To  ke  pmHUked  in  Ike  SHF.  Conference  Troceedingi  (Aug.  1975). 
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Simulation  results  indicate  that  a  bit  rate  reduction  to  1.5  bits/pixel  c»"  be .0b5hiied 

h  .-sst-Sb:  as*s  SS3&M4 

making  the  transform  coding  system  adaptive. 

Predictive  Image  Coding 

The  high  degree  of  correlation  between  a  given  pixel  value  and  its  nearest  neighbors 
allows  an^ image  to  be  efficiently  represented  by  coding  only  the  difference  between  each 
pixel  value  and  its  predicted  value.  The  predicted  pixel  value  is  based  on  previously 
scanned  pixel  values.  In  a  differential  pulse  code  modulation  (DPCM)  system,  the  pre¬ 
dicted  value  of  each  pixel  is  subtracted  from  its  actual  value  and  thLS  difference  is 
quantized  and  transmitted.  Quantization  of  the  data  is  d°^.wlth  a  System 

the  probability  density  of  the  difference  signal.  Thus,  coding  by  use  i of  a  DPCM  system 
requires  a  knowledge  of  the  data  statistics.  The  basic  operation  of  the  DPCM  coder  is  to 
generate  an  uncorrelated  signal  which  is  then  encoded  by  a  memory  less  quantizer 
mission.  At  the  receiver,  the  quantized  difference  signal  is  combined  with  its  predicted 
value  to  form  the  reconstructed  pixel  value.  Basic  DPCM  image  coding  systems  provide  g 
quality  at  about  3  bits/pixel.  Adaptive  systems,  m  which  parameters  of  the  quantizer  and 
predictor  adapt  to  the  image  content,  require  about  2  bits/pixel . 1 51 

Hybrid  Image  Coding 

Analvsis  of  transform  and  predictive  coding  techniques  has  shown  that  both  techniques 
possess  ^tractive  characteristics  and  certain  limitations.  Transform  coding  techniques 
achieve  good  performance  at  low  bit  rates,  show  less  sensitivity  to  data 

ture-to-picture  variations)  and  are  less  vulnerable  to  channel  noise. . 

systems  are  superior  to  transform  techniques  in  terms  of  equipment  complexity,  memory 
requtreme^s?  and  performance  at  high  bit  rates.  Some  limitations  of  predictive  systems 
are  their  sensitivity  to  data  statistics  and  to  channel  error. 

An  intraframe  hybrid  coding  system  that  combines  the  attractive  features  of  {"^trans¬ 
form  and  predictive  coding  systems  has  been  devised. 161  In  this  system,  a  one-dimensional 
transformPis  followed  by  DPCM  linear  predictive  coding  of  the  transform  domain 
cients.  At  the  receiver,  the  transform  coefficients  are  decoded  and  a  replica  of  the 
original  image  is  reconstructed  by  an  inverse  transform. 

Hybrid  Interframe  Image  Coding 

Interframe  coding  of  digital  image  sequences  encompasses  those  techniques  which  make 
use  of  the  high  correlation  that  exists  between  pixel  amplitudes  in  sucq*®®^®,^ra^* 
Intraframe  coding  techniques  that  exploit  spatial  correlations  can,  in  Principle,1 be 
extended  to  include  correlations  in  the  temporal  domain.  Previous  research  in  the  area 
three-dimensional  Fourier  and  Hadamard  transformations  has  indicated  that  bit  rates  can  be 
reduced  by  a  factor  of  five  by  incorporating  correlations  in  the  temporal  direction. 
However ,  three-dimensional  transform  systems  are  unattractive  as  they  use  large  amounts  of 
data  storage  and  require  excessive  computations. 

To  alleviate  the  problems  associated  with  three-dimensional  transform  systems,  new 
hybrid  (two-dimensional  transform) /DPCM  image  coding  systems  have  been  developed. I  1 
These  systems  utilize  both  spatial  and  temporal  correlations  while  greatly  reducing  mem 
ory  sto«gr’nd  computational  requirements.  The  block  diagram  for  a  hybrid  (two-dimen¬ 
sional  transform) /DPCM  system  is  shown  as  figure  1.  In  present  implementations  of  th 
system,  either  a  two-dimensional  cosine  or  Fourier  transformation  is  performed  on  16  x  16 
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Fig.  1.  Hybrid  (two-dimensional  transform) /DPCM  coder. 
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subblocks.  DPCM  linear  predictive  coding  in  the  temporal  domain  is  then  applied  to  the 
transform  coefficients  of  each  subblock.  For  notational  convenience,  the  hybrid  inter¬ 
frame  coders  employing  two-dimensional  Fourier  transforms  will  be  denoted  as  FFD  and  those 
using  two-dimensional  cosine  transforms  as  CCD.  The  FFD  and  CCD  coders  are  adaptive  in 
the  sense  that  statistics  of  the  transform  coefficient  differences  of  each  subblock  are 
computed  prior  to  encoding  the  transform  coefficients  in  the  temporal  direction  by  paral¬ 
lel  DPCM  coders.  At  the  receiver,  the  transmitted  transform  coefficients  are  decoded  and 
a  replica  of  each  frame  is  reconstructed  by  the  appropriate  inverse  two-dimensional  trans¬ 
formation.  These  systems  require  only  a  single  frame  of  storage  and  involve  significantly 
less  memory  and  fewer  computations  than  three-dimensional  transform  coding  techniques. 

Operational  Modes 

At  least  three  operational  modes  have  been  identified  for  the  hybrid  interframe  coding 
systems.  These  operational  modes  depend  on  the  initial  conditions  assumed  for  the  pre¬ 
dictive  coder.  The  initial  conditions  are: 

a.  No  apriori  information  available  at  the  receiver, 

b.  Limited  information  (such  as  mean,  variance  and  temporal  correlations  based  on  a 
statistical  model)  available  at  the  receiver,  and 

c.  First  frame  available  at  the  receiver. 

In  the  no  apriori  information  available  case,  several  frames  are  required  for  the 
hybrid  coder  to  settle.  However,  it  has  been  experimentally  verified  that  in  the  remain¬ 
ing  two  cases,  nearly  stable  coder  performance  is  achieved  within  the  first  4  to  6  frames. 
From  operational  considerations,  the  third  set  of  initial  conditions  is  the  most  realistic 
as  periodic  full  frame  updating  will  be  required  to  eliminate  the  cumulative  effects  due 
to  channel  noise. 

Mathematical  Formulation 

Let  f(x,y)  denote  a  two-dimensional  array  of  intensity  values  on  an  NxN  subblock  of  a 
digital  television  image  of  size  MxM.  Typical  values  for  M  and  N  are  256  and  16,  respec¬ 
tively.  Let  F(u,v)  be  the  two-dimensional  array  obtained  by  taking  the  two-dimensional 
transform  of  f(x,y).  In  the  case  of  the  two-dimensional  discrete  Fourier  transform,  the 
expressions  relating  f(x,y)  and  F(u,v)  are 
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where  N  is  the  size  of  the  square  subblock,  f(x.y)  is  the  image  intensity  function  in  the 
•patial  domain,  x  and  y  are  spatial  coordinates,  F(u,v)  is  the  Fourier  transform,  and  u 
and  v  are  spatial  frequencies. 

For  image  processing  applications,  f(x,y)  is  a  positive  real  function  representing 
brightness  of  the  spatial  sample.  The  two-dimensional  Fourier  transform  of  a  real-valued 
function  has  the  following  conjugate  symmetry  property: 

F* (u,v)  -  F(N-u,N-v)  ,  u,v  »  1,  2,  ....  j  -  1 

The  Fourier  transform  consists  of  2N^  components,  i.e.,  the  real  and  imaginary  or  magni¬ 
tude  and  phase  components  of  each  spatial  frequency.  However,  as  a  result  of  the  con¬ 
jugate  symmetry  properties  mentioned  above,  only  N2  components  are  required  to  completely 
define  the  Fourier  transform. 1*1 

In  the  case  of  the  Fourier  transform,  a  shift  in  the  spatial-domain  variables  results 
in  a  multiplication  of  the  Fourier  transform  of  the  un-shifted  image  by  a  phase  factor. 

If  the  input  image  f(x,y,ti)  is  shifted  by  the  amount  xq  in  the  x-direction  and  y<)  in  the 
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•/-direction  between  times  t!  and  t2,  then  the  Fourier  transform  of  the  shifted  image  is 
given  by 


F(u,v,t2)  «=  F(u,v,tj)  exp  <uxQ  +  vyQ)| 


This  shifting  property  is  expected  to  be  useful  in  detecting  and  compensating  for 
effects  of  motion  between  frames  since  many  types  of  motions  such  as  panned  motion  produce 
significant  changes  in  phase  components  and  small  changes  in  amplitude  components.  Thus, 
compensation  for  platform  motion  may  be  implemented  directly  in  the  array  of  phase  compo¬ 
nents  by  application  of  appropriate  phase  correction  factors. 

The  two-dimensional  Fourier  transform  F(u,v)  of  a  spatial  signal  function  f(x,y)  is 
separable,  i.e.,  it  can  be  computed  as  two  sequential  one-dimensional  transforms  since  the 
Fourier  kernel, 

exp  |i  (ux  +  vy)|  , 

is  separable  and  symmetric.  Thus,  the  basic  one-dimensional  discrete  Fourier  transform 
that  must  be  performed  is 


l  N"1  /  2  i  \ 

F(u>  =  fj  Z  f(x)  exp  ( - N~  ux) 


for  u=0,  1,  .  .  N-l. 

In  the  case  of  the  discrete  Cosine  transforml 3 ] ,  the  one-dimensional  transform  is 


1  N_1 

Flill  '  g  £  f  (x)  cos 


(2x+l)un 

ZN 


for  u=0,  1,  .  .  N-l.  The  form  of  eq.  (5)  differs  from  that  of  reference  13]  only  by  a 
normalization  constant.  The  cosine  transform  is  also  separable  and  a  two-dimensional  dis¬ 
crete  cosine  transform  of  an  NxN  subblock  results  in  N  real  coefficients. 

Quantization 

Experimental  evidence  derived  from  transmission  of  a  typical  "head  and  shoulders”  pic¬ 
ture  telephone  scene  has  shown  that  the  frame  difference  signal  has  a  probability  density 
closely  approximated  by  a  double  sided  exponential  function. 1 101  The  optimum  minimum 
mean  square  error  quantizer  for  this  distribution  has  been  found  to  be  a  uniform  quantizer 
combined  with  a  companding  of  the  frame  difference  signal.  I11' 

Since  the  variances  of  the  transform  domain  coefficient  different  'S  are  different,  it 
is  necessary  to  use  different  quantizer  parameters  for  each  one.  Each  coefficient  differ¬ 
ence  signal  is  allocated  a  number  of  bits  proportional  to  its  estimated  variance  in 
accordance  with  an  optimum  bit  assignment  algorithm. 

Fidelity  Criteria 

In  figure  1,  differences  between  input  signal  f(x,y,t)  and  output  signal  f(x,y,t) 
are  due  to  two  sources:  quantization  errors  and  errors  due  to  channel  noise.  To  evalu¬ 
ate  coding  efficiency  of  the  hybrid  encoders,  two  objective  criteria  were  used.  The  first 
criterion,  NMSE,  is  a  measure  of  the  mean  square  error  between  f(x,y,t)  and  f(x,y,t) 
averaged  over  an  entire  frame  of  size  MxM.  Normalization  is  achieved  by  dividing  the  mean 
square  error  by  the  mean  signal  energy  within  the  frame. 

M-l  M-l  ,  „  I, 

Z  Z_  -  f ix,y,t)| i 


f (x,y, t) 


The  second  criterion,  SNR,  measures  the  ratio  of  peak-to-peak  signal  to  RMS  noise: 
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Figures  2  and  3  are  graphs  illustrating  the  coding  efficiency  of  the  hybrid  FFD  and  CCD 
coders  at  various  bit  rates  in  the  interval  0.1  to  1.0  bits/pixel/frame.  To  perform  this 
series  of  experiments,  a  256  x  256  resolution  data  base  consisting  of  16  consecutive 
frames  of  a  24  frames  per  second  (fps)  motion  picture  was  digitized.  Initial  conditions 
assumed  were  that  the  first  frame  was  available  at  the  receiver. 
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Fig.  2.  Four ier/Fourier/DPCM  coder 
at  various  bit  rates. 


Fig.  3. 


Cosine/cosine/DPCM  coder 
at  various  bit  rates. 


Photographs  of  frame  number  16  after  coding  by  the  FFD  and  CCD  coders  at  average  pixel 
bit  rates  of  1.0,  0.5,  0.25,  and  0.1  are  shown  as  figures  4  and  5.  The  results  shown  in 
figure  4  for  the  FFD  coder  were  obtained  by  coding  the  real  and  imaginary  components  of 
the  Fourier  coefficients  by  assigning  half  of  the  available  bits  to  each  component. 


m-,%  m-% 

1.0  bits/pixel/frame  0.5  bits/pixel/frame  0.25  bits/pixel/frame  0.1  bits/pixel/f rame 

Fig.  4.  FFD  coder  for  frame  16. 


1.0  bits/pixel/frame  0.5  bits/pixel/frame  0.25  bits/pixel/frame  0.1  bits/pixel/frame 

Fig.  5.  CCD  coder  for  frame  16. 
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Noise  Immunit 


Pei formance  of  the  FFD  and  CCD  hybrid  interframe  coders  was  investigated  in  the  pres¬ 
ence  of  channel  noise.  In  order  to  study  the  effect  of  channel  noise,  a  binary  symmetric 
channel  was  simulated.  The  channel  is  assumed  to  operate  on  each  binary  digit  independ¬ 
ently,  changing  each  digit  from  0  to  1  or  from  1  to  0  with  probability  PCe  and  leaving  the 
digit  unchanged  with  probability  1-Pge.  At  the  receiver,  the  encoded  picture  is  recon¬ 
structed  from  the  string  of  binary  digits,  including  errors,  transmitted  across  the  channel. 

Dearadations  due  to  channel  noise  probabilities.  Pee'  zero,  10  ^  and  10  ^  for  the 
FFD  and  CCD  coders  at  average  bit  rates  of  1.0  and  0.25  bi ts/pixel/f rame  are  shown  in 
figures  6  through  9.  The  generally  monotonical ly  increasing  character  cf  these  curves 
illustrates  the  fact  that  once  an  error  has  occurred,  it  tends  to  propagate  in  the  tem¬ 
poral  direction  until  corrected  by  a  frame  refresh. 


Fig.  6.  Effects  of  channel  noise  for  Fig.  7.  Effects  of  channel  noise  for 

Four ier /Four ier/D PCM  coder  at  Fourier/Fourier/DPCM  coder  at 

an  average  1.0  bits/pixel/frame.  an  average  0.25  bits/pixel/frame. 


Fig.  8.  Effects  of  channel  noise  for  Fig.  9.  Effects  of  channel  noise  for 

cosi ne/cosine/DPCM  coder  at  an  cosine/cosine/DPCM  coder  at  an 

average  1.0  bits/pixel/frame.  average  0.25  bits/pixel/frame. 


Photographs  corresponding  to  average  bit  rates  of  1.0  and  0.25  bits/pixel/frame  for  the 
FFD  and  CCD  coders  with  channel  error  probabilities  of  10"  3  and  10"^  are  shown  in  figures 
10  and  11. 
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Bit  Transfer  Rate 


In  keeping  with  the  previously  mentioned 
transmitted  while  retaining  image  fidelity, 
certain  bit  transfer  rates  (BTR)  across  the 
product  of  average  pixel  bit  rate  per  frame 


objective  of  minimizing  the  number  of  bits 
a  series  of  experiments  was  performed  in  which 
channel  were  fixed.  The  BTR  is  defined  as  the 
and  frame  rate  and  has  units  of  bits/pixel/sec 


BTR  =  (bits/pixel/frame)  *  (frames/sec)  (8) 

The  available  16  frame  test  data  base  was  extracted  from  a  24  fps  motion  picture 
sequence.  By  employing  frame  skipping  techniques,  temporal  subsampling  was  used  to  simu¬ 
late  short  12,  8  and  6  fps  sequences  from  the  16  frame  test  data  base. 


Average  bit  rates  in  the  interval  0.083  to  1.333  bits/pixel/frame  were  used  in  conjunc¬ 
tion  with  the  four  frame  rates  mentioned  above  to  perform  simulations  with  BTR  values  of 
8,  6,  4  and  2  bits/pixel/sec.  The  results  of  these  experiments  are  shown  as  figures  12 
through  15,  respectively.  For  all  cases  examined,  the  graphs  show  that  reduced  frame  rates 
produce  smaller  NMSE  values  for  the  individual  frames  coded.  This  indicates  that  reduc¬ 
tions  experienced  in  f rame-to-f rame  correlations  due  to  temporal  subsampling  are  completely 
compensated  for  by  th  •  increased  number  of  bits  available  for  coding.  However,  subject¬ 
ively,  reduced  frame  rates  tend  to  result  in  jerky  subject  motion.  This  is  most  apparent 
for  rapidly  moving  objects  in  the  field  of  view  and  is  of  lesser  consequence  for  slowly 
changing  scenes. 
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Fig.  12.  Cosine/cosine/DPCM  coder  at 
bit  transfer  rate  of  8  bits/ 
pixel/sec. 


Fig.  13.  Cosine/cosine/DPCM  coder  at 
bit  transfer  rate  of  6  bits/ 
pixel/sec. 
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Fig.  14.  Cosine/cosine/DPCM  coder  at 
bit  transfer  rate  of  4  bits/ 
pixel/sec. 


rig.  15.  Cosine/cosine/DPCM  coder  at 
bit  transfer  rate  of  2  bits/ 
pixel/sec . 


Conclusions 

Based  on  theoretical  and  experimental  results  obtained  to  date,  two  main  conclusions 
have  been  reached.  The  first  is  that  exploitation  of  temporal  correlations  in  addition  to 
spatial  correlations  has  been  demonstrated  to  be  a  viable  technique  for  coding  sequences 
of  digital  images.  This  fact  is  demonstrated  by  a  comparison  of  the  average  bit  rates 
required  for  the  interfr.ime  cosine/cosine/DPCM  and  the  existing  intraframe  cosine/DPCM 
coders  to  achieve  the  same  level  of  NMSE  performance.  The  sixteenth  frame  of  the  test 
data  base  was  chosen  for  comparison  and  was  coded  at  an  average  0.25  bits/pixel  by  the 
interframe  cosine/cosine/DPCM  coder.  When  using  the  intrafrane  cosine/DPCM  coder,  it  was 
necessary  to  code  this  frame  at  a  bit  rate  of  more  than  2  bits/pixel  to  achieve  the  same 
NMSE. 

The  second  conclusion  is  that  the  performance  of  the  hybrid  interframe  coders  investi¬ 
gated  are  heavily  data  dependent.  In  the  case  of  the  16  frame  head  and  shoulders  test 
data  base,  good  coding  performance  was  achieved  since  subject  movement  was  restricted  to  a 
relatively  small  portion  of  the  image.  However,  coding  performance  with  a  different  aerial 
data  base  was  degraded  from  the  previous  case  due  principally  to  platform  motion  which 
caused  f rame-to-f rame  pixel  amplitude  variations  across  the  entire  image.  Since  the  per¬ 
formance  of  the  hybrid  interframe  coders  is  dependent  on  temporal  correlation,  a  reduced 
level  of  performance  is  to  be  anticipated  for  image  sequences  distorted  by  motion. 


58 


References 


1.  P.  A.  Wintz,  ''Transform  Picture  Codinq,"  Proceedinqs  IEEE,  Vol.  60,  No.  7. 
pp.  809-820,  July  1972. 

2.  W.  K.  Pratt,  L.  R.  Welch  and  W.  Chen,  "Slant  Transforms  for  Image  Coding," 

IEEE  Transactions  on  Communications.  Vol.  COM-22,  No.  8,  pp.  1075-1093,  August,  1974. 

3.  N.  Ahmed,  T.  Natarajan,  K.  R.  Rao,  "Discrete  Cosine  Transform,"  IEEE  Trans,  on 
Computers,  Vol.  C-23,  No.  1,  pp.  90-93,  January  1974. 

4.  W.  K.  Pratt,  "Spatial  Transform  Coding  of  Color  Images,"  IEEE  Transactions  on 
Communication  Technology,  Vol.  COM-19,  No.  6,  pp.  980-992,  December,  1971. 

5.  A.  Habibi,  “Comparison  of  nth  -order  DPCM  Encoder  with  Linear  Transformation  and 
Block  Quantization  Techniques,"  IEEE  Transactions  on  Communication  Technology,  Vol.  COM- 
19,  No.  6,  Part  I,  pp.  948-956,  December  19^1. 

6.  A.  Habibi,  "Hybrid  Coding  of  Pictorial  Data,"  IEEE  Transactions  on  Communications, 
Vol.  COM-22,  No.  55,  pp.  614-624,  May  1974. 

7.  A.  G.  Tescher,  "The  Role  of  Phase  in  Adaptive  Image  Coding,"  Ph.D.  Thesis, 
University  of  Southern  California,  Electrical  Engineering  Department,  January  1974.  Pub¬ 
lished  as  Report  510,  University  of  Southern  California,  Image  Processing  Institute. 

8.  J.  A.  Roese,  W.  K.  Pratt,  G.  S.  Robinson  and  A.  Habibi,  "Interframe  Transform 
Coding  and  Predictive  Coding  Methods,"  Proceedinqs  of  1975  International  Conference  on 
Communications  (ICC  75),  Vol.  II,  Paper  i3,  ppT  17-21,  June  16-18,  1975T 

9.  G.  S.  Robinson,  "Orthogonal  Transform  Feasibility  Study,"  NASA  Final  Report  NASA- 
CR-115314,  N72-13143  (176  pp. )  (submitted  by  COMSAT  Labs,  to  NASA  Manned  Spacecraft 
Center,  Houston,  Texas)  November  1971. 

10.  A.  J.  Seyler,  "Probability  distributions  of  television  frame  differences,"  Proc. 

I  REE,  Australia,  pp.  355-366,  November  1965. 

11.  B.  Smith,  "Instantaneous  companding  of  quantized  signals,"  Bell  System  Technical 

Journal,  Vol.  36,  pp.  653-709,  May  1957.  "  - 

Acknowledgement 

This  research  was  supported  by  the  Advanced  Research  Projects  Agency  of  the  Department 
of  Defense  and  monitored  by  the  Air  Force  Eastern  Test  Range  under  Contract  R08606-72- 
0008  and  by  the  Naval  Undersea  Center,  San  Diego,  California,  under  Contract  N00123-75-C- 


59 


1 


APPENDIX  C 

INTERFRAME  TRANSFORM  CODING 
AND  PREDICTIVE  CODING  METHODS 


INTERFRAME  TRANSFORM  CODING  AND  PREDICTIVE  COOING  METHODS  * 

John  A.  Roese  Ali  Habibi 

Naval  Undersea  Center  TRW  Systems  Group 

San  Diego,  California  92 1 32  Redonao  Beach.  California  90278 

William  K.  Pratt  Guner  S.  Robinson 

Image  Processing  Institute 
University  ol  Southern  California 
Los  Angeles.  California  90007 


ADSTRAC-T 

Transform,  differential  pulse  code  modulation 
(DPCM),  and  hybrid  transform/DPCM  image  coding 
methods  are  applied  to  coding  successive  frames  of 
digital  images.  These  coding  techniques  are  design¬ 
ed  to  exploit  the  inherent  spa  al  and  temporal  cor¬ 
relations  of  image  sequences.  In  transform  image 
coding,  subsections  of  images  are  transformed  into 
arrays  of  nearly  uncorrelated  coefficients  by  the  use 
of  orthogonal  transformations.  Bit  compression  is 
achieved  by  assigning  the  available  bits  in  proportion 
to  the  energy  of  the  transform  coefficients.  In 
DPCM  image  coding  systems,  the  value  of  the  image 
sample  is  predicted  and  the  difference  between  the 
actual  and  the  predicted  value  is  quantized  and  trans¬ 
mitted.  Hybrid  transform/DPCM  coding  implemen¬ 
tations  combine  the  relatively  superior  low  bit  rate 
performance  and  channel  noise  immunity  of  trans¬ 
form  methods  with  the  minimal  storage  requirements 
of  DPCM  encoders.  Interframe  coding  implemen¬ 
tations  studied  include:  three-dimensional  Cosine 
transforms,  hybrid  two-dimensional  spatial  trans¬ 
forms  with  DPCM  in  the  temporal  direction,  e.  g.  , 
Cosine-Cosine/ DPCM  and  Fourier-Fourier /DPCM, 
and  three-dimensional  DPCM.  Results  are  sum¬ 
marized  by  evaluating  the  various  implementations 
in  terms  of  performance  and  complexity  criteria. 

INTRODUCTION 

During  the  past  few  years  several  intraframe 
digital  image  coding  systems  based  upon  transform 
and  linear  predictive  coding  concepts  have  been  de¬ 
veloped  1,2  .  Most  of  these  systems  have  achieved 
bit  rate  reductions  by  removal  of  statistical  redun¬ 
dancies  within  an  image  frame  combined  with  the  de¬ 
letion  of  that  part  of  the  spatial  image  representation 
least  critical  to  the  human  observer.  It  is  known 
that  there  is  considerable  temporal  redundancy  be¬ 
tween  frames  of  real  time  imagery  systems;  also 
there  are  psychophysical  limitations  to  temporal 
image  perception.  Exploitation  of  either  property 
in  the  past  has  been  difficult  to  achieve  because  of 
implementation  problems.  However,  there  have 
been  recent  technological  advances  in  signal  proces¬ 
sing  circuitry  which  hold  promise  for  the  practical 
implementation  of  digital  real  time  television  image 
coders.  Several  such  systems  are  described  in 
this  paper  along  with  an  analysis  of  their  perfor¬ 
mance. 

INTRA  FRAME  TRANSFORM  AND  PREDICTIVE 

CODING 

Transform  coding  and  linear  predictive  coding 


mTo  be  published  lit  ICC  Conference  Proceedings  (Jen.  1975). 


are,  perhaps,  the  most  widely  employed  techniques 
for  intraframe  image  coding.  Operational  descript¬ 
ions  of  these  coding  methods  are  given  below. 

Transform  Coding:  In  transform  coding  systems, 
a  one  or  two-dimensional  mathematical  transform  of 
an  image  line  segment  or  block  is  performed.  The 
resulting  transform  coefficients  are  then  quantized 
and  coded.  A  bit  rate  reduction  is  possible  because 
the  distribution  of  energy  in  the  transform  domain 
permits  more  efficient  quantization  and  coding. 

After  decoding,  an  inverse  transform  is  taken  to 
obtain  a  replica  of  the  image  at  the  receiver.  Trans¬ 
forms  that  have  proven  useful  include  the  Fourier, 
Hadamard,  Slant,  Karhunen- Loeve,  and  Cosine  trans¬ 
forms  [3-5\  Transform  coding  techniques  have 
been  explored  extensively  both  theoretically  and  by 
simulation.  It  has  been  shown  that  a  significant  bit 
rate  reduction  can  be  achieved  in  many  applications 
with  minimal  image  degradation.  Simulation  results 
indicate  that  a  bit  rate  reduction  to  1.  5  bits /pixel  can 
be  obtained  for  monochrome  image  transform  coding 
in  16  x  16  pixel  blocks,  while  color  images  require 
about  2.  Obits  pixel  ‘6  .  The  bit  rate  can  be  re¬ 
duced  further  by  making  the  transform  coding 
system  adaptive. 

Predictive  Coding:  The  high  degree  of  corre¬ 
lation  between  a  given  pixel  value  and  its  nearest 
neighbors  allows  an  image  to  be  efficiently  repre¬ 
sented  by  coding  only  the  difference  between  each 
pixel  value  and  its  predicted  value.  The  predicted 
pixel  value  is  based  on  previously  scanned  pixel 
values.  In  a  differential  pulse  code  modulation 
(DPCM)  system,  the  predicted  value  of  each  pixel  is 
subtracted  from  its  actual  value  and  this  difference 
is  quantized  and  transmitted.  At  the  receiver,  the 
quantized  difference  signal  is  combined  with  its  pre¬ 
dicted  value  to  form  the  reconstructed  pixel  value. 
Basic  DPCM  image  coding  systems  provide  good 
quality  at  about  3  bits /pixel.  Adaptive  systems,  in 
which  parameters  of  the  quantizer  and  predictor 
adapt  to  the  image  content,  require  about  2  bits /pixel. 

INTERFRAME  TRANSFORM  CODING 

In  interframe  transform  coding,  a  three-dimen¬ 
sional  unitary  transform  is  performed  on  the  data. 

Let  fix,  y,  t)  denote  a  three-dimensional  array  of 
amplitude  values  for  each  frame  of  a  digital  image. 
Also,  let  F(u,v,  w'  be  the  three-dimensional  array 
obtained  by  taking  the  three-dimensional  transform 
in  the  (x#y,t)  domain.  If  the  size  of  the  three- 
dimensional  array  is  NjxN2*Nj  ,  then  such  a  trans¬ 
form  coder  can  be  described  in  a  general  form  as 
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F(u,v,  w)  =  EEE  f<x,y,t)0lu,v,  w,  x.y.t) 
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where  x  and  y  are  spatial  coordinates,  t  is  the  time 
coordinate,  u,  v,  w  are  the  transform  domain  coordi¬ 
nates,  and  0lu.v,w,  x,y,t)  represents  a  set  of  three- 
dimensional  basis  vectors.  For  the  three-dimen¬ 
sional  discrete  Fourier  transform,  eq.  11)  has  the 


form. 
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f(x,y.t)  =  5353  ^Flu.v.w). 
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exp 

Since  the  Fourier  kernel 


HI?*,*,)] 


exp 


is  separable,  the  three-dimensional  transform  can 
be  computed  as  a  one-dimensional  transform  in  the 
temporal  direction  followed  by  a  two-dimensional 
transform  in  the  spatial  domain.  The  two-dimen¬ 
sional  spatial  transform  car*  in  turn,  be  computed 
as  a  one-dimensional  transform  along  the  rows 
followed  by  a  one-dimensional  transform  along  the 
columns.  Thus,  the  basic  one-dimensional  trans¬ 
form  that  must  be  performed  i  s 


F(u)  = 


-L£f,x)  **p(-^-ux) 

x-  0 


The  form  of  eq.  (4)  differs  from  that  of  reference  [71 
only  by  a  normalization  constant.  An  example  of 
three-dimensional  transform  coding  of  a  16  frame  data 
base  is  given  in  figure  1.  This  figure  illustrates  the 
decoded  images  for  frames  number  1  and  16.  In  this 
example,  a  three-dimensional  Cosine  transform  wis 
performed  on  cubic  bloc  ts  of  size  16x16x16  on  16 
frames  of  size  Z56xZ56.  The  average  bit  rate  used 
was  1.  0  bits/pixel. 


(Z) 


for  u=0, 1, .  .  . ,  N-l.  In  the  case  of  the  discrete 
Cosine  transform  [7  1,  the  one -dimensional  trans 
form  ia 

ri*” 


for  u=  0, 1,  •  •  •  »  N-l# 


F  rame  1 


F  rame  16 


Figure  1 


INTER  FRAME  DPCM  CODING 

In  DPCM  systems,  linear  prediction  is  used  to 
generate  a  differential  signal  which  is  quantized  and 
transmitted.  At  the  receiver  a  similar  predictor 
uses  some  previously  transmitted  values  of  the  quan- 
tized  differential  signal  to  obtain  a  facsimile  of  the 
transmitted  signal. 

Prediction  of  a  picture  element  is  performed  by 
using  a  set  of  previously  scanned  picture  elements 


£ 


(5) 


where  f  is  the  set  of  picture  elements  with  zero 
mean  and  variance  -2,and  n  is  the  order  of  the  pred.c- 
tor.  The  predictor  parameters  a;  are  the  solutions 
of  n  algebraic  equations: 


R  =>  a  R. .  i-  1, 2, 
oi  Z-s  J  0 


(6) 


(3) 


where  Rij  are  the  correlations  between  picture  ele- 
menta: 


(4) 


V 


The  variance  of  the  differential  signal 


(7) 
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INTERFRAME  HYBRID  CODING 


«•  5  -  S 


is  related  to  the  variance  ~  of  picture  elements  by 


n 

i  1 

In  a  three-dimensional  DPCM  interframe  coder, 
the  value  of  the  picture  element  s0  is  estimated  using 
previously  scanned  picture  elements  as  shown  in 
figure  2  in  which  sC)  is  the  picture  element  to  be 
predicted.  The  elements  sj,  S2,  S4,  and  S5  are  the 
previously  scanned  picture  elements  in  the  present 
frame,  and  S3,  Sy,. . .  ,  Sj2>  SJ3  are  the  closest 
samples  on  the  previous  frame. 


Figure  2 

Figure  3  shows  frames  number  2  and  16  of  the 
16  frame  data  base  after  coding  by  a  three-dimen¬ 
sional  third  order  predictive  coder  at  2.  0  bits/pixel 
The  picture  elements  used  for  predicting  s0  are  Sj 
and  S2  on  the  present  frame  and  S3  on  the  previous 
frame. 


Frame  2  Frame  16 

Figure  3 


Three-dimensional  hybrid  encoders  have  been 
investigated  which  use  two-dimensional  transforma¬ 
tions  in  the  spatial  domain  cascaded  with  a  DPCM 
encoder  in  the  temporal  domain.  Due  to  DPCM  en¬ 
coding  in  the  temporal  direction,  these  encoders  re¬ 
quire  storage  of  only  the  single  previous  frame  of 
data.  Simulation  results  indicate  that  the  hybrid 
three-dimensional  encoders  perform  better  than  the 
corresponding  three-dimensional  transform  encoders. 
This  observation  is  confirmed  by  the  comparative- 
performance  of  the  two-dimensional  hybrid  and 
transform  encoders  previously  reported  in  the  lit  — 
erature  8  ,  The  following  discusses  two  different 
three-dimensional  hybrid  coding  algorithms  which 
have  been  investigated. 

Two-Dimensional  Cosine  Transform/ DPCM 
Coder:  This  hybrid  coder  exploits  spatial  image 
correlations  by  taking  a  two-dimensional  discrete 
Cosine  transform  and  temporal  correlation  by  use  of 
a  DPCM  coder.  Theoretical  studies  indicate  that 
this  hybrid  interframe  coder  possesses  the  attrac¬ 
tive  features  of  the  hybrid  intraframe  coder.  It  is 
anticipated  that  this  system  will  reduce  the  number 
of  bits  needed  for  reconstruction  of  television  images 
at  the  receiver  by  a  factor  of  five  over  the  two- 
dimensional  hybrid  coder.  Figure  4  shows  the  results 
for  frames  1  and  16  after  applying  the  hybrid  Cosine- 
Cosine/DPCM  encoder.  The  average  bit  rate  used 
was  1.0  bits/pixel. 


F  ratne  1  F  rame  16 

Figure  4 


Two-Dimensional  Fourier  T ransform/ DPCM 
Coder:  An  alternative  approach  to  the  three-dimen¬ 
sional  hybrid  coder  uses  the  two-dimensional 
Fourier  transform  in  place  of  the  two-dimensional 
Cosine  transform.  In  this  system,  the  two-dimen¬ 
sional  Fourier  transform  coefficients  of  each  frame 
are  coded  and  transmitted.  Figure  5  illustrates  the 
results  for  frames  1  and  16  due  to  Fourier- Fourier  / 
DPCM  coding  of  the  16  frame  data  base.  These 
results  were  obtained  by  coding  the  real  and  imagi¬ 
nary  components  of  the  Fourier  coefficients  at  an 
average  bit  rate  of  0.  5  bits  for  each  component. 
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Frame  1  Frame  1 6 

Figure  ? 

The  two-dimensional  Fourier  transform  of  the 
kth  frame,  Fnfu,  v),  can  be  expressed  as 


1*-  (U.  v 

F,  iu,  vl  A,  iu,  vie*  k 
k  k 


where  Afciu.v)  and  v)  are  the  amplitude  and  phase 

planes  of  the  kth  frame.  Many  types  of  motions  such 
as  panned  motion  produces  significant  changes  in  the 
phase  plane  and  small  changes  in  the  amplitude  plane. 
An  attractive  feature  of  amplitude-phase  coding  is 
that  motion  compensation  may  be  implemented 
directly  in  the  phase  plane  by  application  of  appro¬ 
priate  phase  correction  factors. 

SY ST EMS  ANALYSIS 


In  this  section,  direct  comparisons  are  made 
between  the  four  previously  discussed  interframe 
encoders.  The  comparison  criteria  used  are:  coding 
efficiency,  storage  requirements,  noise  immunity 
and  implementation  complexity. 

The  results  shown  as  figures  1,3,  4,  and  5  il¬ 
lustrate  the  coding  performance  of  the  four  encoders. 
For  MSE  comparison  purposes,  the  encoders  were 
run  at  an  average  bit  rate  of  1.0  bits  /pixel.  This 
bit  rate  was  chosen  as  it  is  the  lowest  bit  rate  at 
which  the  three-dimensional  DPCM  encoder  can 
operate.  Additional  experiments  have  shown  that  the 
hybrid  transform/DPCM  encoders  can  be  successfully 
operated  at  even  lower  bit  rates.  For  each  encoder, 
MSE  calculations  were  made  for  each  frame  coded. 
The  MSE  values  were  normalized  relative  to  total 
image  energy  of  each  frame.  Comparison  of  the 
normalized  MSE  values  indicate  that  the  hybrid 
Cosine- Cosine/DPCM  and  P  ourier- Fourier /DPCM 
encoder  implementations  were  superior  in  terms  of 
MSE  coding  efficiency.  This  conclusion  is  supported 
by  subjective  comparison  of  the  coded  frames  illus¬ 
trated  above. 

A  significant  disadvantage  of  all  interframe  cod¬ 
ing  systems  is  the  requirement  for  storage  of  pre¬ 
viously  scanned  data  frames.  Of  the  interframe 
systems  considered  in  this  paper,  the  three-dimen¬ 
sional  Cosine  implementation  is  the  least  attractive 
in  terms  of  required  storage  as  several  previous 
frames  must  be  retained.  Even  if  the  number  of 


fraiiU'S  stored  is  const  mined  ti<  a  ln'.v  number,  such 
•.a  iu-ir,  tin-  total  memory  requirements  are  still 
excessive.  This*  is  evident  when  the  three- dimen¬ 
sional  Cosine  encoder  is  compared  to  the  three  other 
implementations  which  use  first-order  predictive 
DPCM  coding  in  the  temporal  domain.  Under  the 
assumption  of  first-order  temporal  predictors,  only 
the  single  previous  frame  needs  to  be  stored.  A 
special  case  occurs  for  the  three-dimensional  DPCM 
implementation  where,  in  addition  to  the  previous 
frame,  it  is  also  necessary  to  retain  the  previously 
scanned  line  of  the  current  frame. 

Immunity  to  channel  noise  varies  widely  for  the 
interframe  coder  implementations  considered.  The 
least  sensitive  is  the  three-dimensional  transform 
encoder,  whereas  the  DPCM  encoder  is  most  vainer- 
able  to  channel  noise  due  to  its  transmission  of 
simple  pixel  amplitude  differences  at  a  low  bit  rate. 
The  hybrid  transform/predictive  encoders  transmit 
differences  in  transform  coefficients  instead  of  pixel 
amplitudes  and  are  less  sensitive  to  channel  noise 
than  strictly  predictive  encoders.  An  aspect  of 
predictive  interframe  coders  which  has  yet  to  be 
fully  investigated  is  their  amenability  to  use  with 
error  detection  and  correction  algorithms. 

The  implementation  complexity  criterion  is  a 
coarse  measure  of  weight,  size  and  power  require¬ 
ments  for  each  encoder.  The  inherently  simple  oper¬ 
ation  of  DPCM  coders  combined  with  essemiallv  a 
single  frame  of  storage  favors  the  DPCM  interframe 
coder.  Conversely,  the  multiple  frame  storage 
requirements  of  the  three-dimensional  transform 
coders  severely  limit  their  usefulness.  The  two 
hybrid  transform/predictive  encoders  appear  to  be 
the  best  compromise  as  they  combine  single  frame 
storage  requirements  with  the  simplicity  of  DPCM  in 
the  temporal  domain. 

A  summary  of  the  results  of  the  systems  analysis 
for  the  interframe  encoders  is  contained  in  Table  1. 
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APPENDIX  D 


THE  PRIME  COSINE  TRANSFORM 


'll  If  PRIME  COSINE  TRANSFORM 
J.  Speiser 

NAVAL  UNDERSEA  CENTER 
SAN  DIEGO,  CA  92132 


This  note  discusses  a  high  speed  implementation  of  the  odd  discrete  cosine 
transform  (ODCT)  which  eliminates  the  multipliers  required  in  earlier  imple¬ 
mentations  [1]  based  on  the  chirp-Z  transform.  The  discrete  cosine  transform 
is  useful  for  television  data  compression  since  its  basis  vectors  closely  approx¬ 
imate  those  of  the  optimum  Karhunen-Loeve  transform  for  exponentially  corre¬ 
lated  data  [2]. 

The  ODCT  is  defined  as  the  first  N  Fourier  coefficients  of  the  length  2N-1 
even  extension  of  the  data,  assuming  that  the  data  consists  of  N  real  values. 

This  is  shown  in  eqn  (1). 


G 


k 


N-l 

Z 

n=  -  (N-l) 


-i2irkn/(2N-l) 


for  k=0,l 


.N-l 


(1) 


where  g_n  =  gn 

In  order  to  be  able  to  use  a  variant  of  Rader's  Prime  Transform  algorithm 
[3]  we  assume  that  P  =  2N-1  is  a  prime.  Note  that  the  data  block  length,  N, 
need  not  be  a  prime,  as  shown  in  the  Table. 

It  has  been  previously  shown  that  if  P  is  a  prime,  the  discrete  Fourier 
transform  (DCT)  of  length  P  can  be  implemented  using  a  circular  convolution  of 
length  P-1  together  with  two  analog  permuter  memories  of  length  P  [4],  It  will 
be  shown  here  that  the  symmetry  of  the  extended  data  in  equation  (1)  permits 
the  size  of  the  circular  convolution  and  the  permuter  memories  to  be  reduced  by 
a  factor  of  two. 


For  each  prime  P,  there  is  an  integer  R,  called  a  primitive  root  of  P, 
such  that  the  residues  of  R,  R^,  ...  rP-1  are  all  distinct  modulo  P,  and  include 
every  nonzero  residue  modulo  P  [5],  Therefore,  for  each  integer  n  not  congru¬ 
ent  to  zero  mod  P,  n  can  be  represented  uniquely  as  a  power  of  R  modulo  P,  say 
n  =  Rn  (mod  P).  The  integer  n'  is  called  the  index  of  n  (mod  P)  with  respect 
to  the  primitive  root  R.  In  effect,  R  plays  the  role  of  the  base  of  a  system 
of  logarithm  in  modulo  P  arithmetic,  and  n'  is  the  logarithm  of  n.  This  repre¬ 
sentation  is  useful  because  it  allows  us  to  replace  multiplication  by  addition 
in  the  exponent  of  the  DFT,  and  thus  reduce  the  DFT  to  a  circular  correlation, 
as  shown  in  equation  (2). 
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Since  zero  does  not  have  an  index  (logarithm)  with  respect  to  the  primi¬ 
tive  root,  the  zero  frequency  point  in  the  transform  must  be  computed  sepa¬ 
rately,  as  shown  in  equation  (3). 
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(3) 


The  interpretation  of  equation  (2)  is  that  the  DFT  coefficients  in  per¬ 
muted  order  are  obtained  by  adding  gQ  to  the  circular  crosscorrelation  of  a 
permuted  sinusoid  with  a  permuted  version  of  the  data  points  excluding  gp. 

For  the  special  case  of  the  DCT,  the  symmetry  of  the  extended  data  allows  us 
to  replace  the  complex  exponential  by  a  cosine,  as  shown  in  equation  (4). 

N_1  v'+n' 

CRk'  =  gQ  +  £  gRn'  cos  (2ttRk  n  /  (2N-1) )  (4) 

n'=  -(N-l) 

n'^  0 


It  will  now  be  shown  that  the  permuted  data  and  the  permuted  cosine  have 
periodicity  N-l,  so  that  the  circular  correlation  of  length  P-1  =  2N-2  can  be 
reduced  to  a  circular  correlation  of  length  N-l.  First,  note  that  both  the 
extended  data  and  the  cosine  function  are  even.  Let  h  be  any  even  sequence. 
It  will  be  shown  that  hRS  has  period  N-l,  where  the  subscript  Rs  is  reduced 
modulo  P  =  2N-1. 

It  is  well  known  in  number  theory  that  =  -1  (mod  P)  [5].  In  our 

case,  (P-l)/2  =  (2N-2)/2  =  N-l.  Therefore  hRS  +  (N-l)  =  hRs  RN-1  =  h_RS  =  hRs 

Using  this  periodicity  property  applied  to  equation  (4)  lets  us  write  the 
ODCT  in  shorter  form,  as  shown  in  equation  (5) . 
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(5A) 
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(5B) 


The  circular  correlation  required  for  equation  5A  may  be  implemented  by 
any  of  the  alternative  methods  for  implementing  a  circular  correlation.  The 
most  straightforward  would  be  to  use  a  transversal  filter  of  length  2(N-1)-1  = 
2N-3,  with  tap  weights  of  cos (2TOS/ (2N-1) ) ,  for  s  =  N-l , . . . 1 ,0, 1 , 2 . . .N-l . 

The  architecture  of  the  transform  is  shown  in  Figure  1,  and  is  virtually 
identical  to  that  of  the  prime  Fourier  transform,  except  that  the  prime  cosine 
transform  need  only  permute  real  data  and  filter  the  permuted  real  data  with  a 
filter  having  real  weights.  The  prime  cosine  transform  is  thus  considerably 
simpler  to  implement  than  a  prime  Fourier  transform  of  the  same  block  length. 
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A  suitable  analog  permuter  memory  has  been  developed  for  the  Naval  Undersea 
Center  as  a  minor  modification  of  a  commercially  available  serial  analog  memory*. 
The  commercially  available  serial  access  memory  stores  analog  samples  as  charges 
in  an  array  of  MOS  capacitors  under  the  control  of  read  in  and  read  out  shift 
registers.  The  permuter  memory  shown  in  figure  2  differs  only  in  the  fact  that 
one  of  the  shift  registers  has  been  replaced  by  a  binary  decoder,  thus  allow¬ 
ing  the  data  to  be  reordered  by  an  external  control  signal. 


*Reticon  Corp.  Sunnyvale,  CA. 
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p  (prime) 


N  =  (P+l)/2  (data  block  length) 


2N-3  =  P-2  =  filter  length 


31 

16 

29 

61 

31 

59 

127 

64 

125 

2S1 

126 

249 

2S7 

129 

255 

Table  1 


Selected  Primes  and  the  corresponding  OCDT  lengths 
and  filter  lengths.  (The  filter  lengths  shown 
assume  that  the  data  is  not  recirculated 
or  reread  into  the  filter.) 


SAMPLE  AND  HOLD 


n  (HOLD  FOR  N  time  units) 


INTEGRATOR 
(INTEGRATE  FOR  N-l  time  units)! 


C  vtATi"Gi 


PERMUTER 

— 

FILTER 

»J  PERMUTER 

i _ , _ _ 

(the  switch  is  in  the  up  position  for  the  first  data  sample,  and  is  down  for 
the  remaining  N-l  samples) 

Figure  1.  Prime  ODCT  architecture. 
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FIGURE  2 

ANALOG  PERMUTATION  MEMORY 
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PRIME  TRANSFORM  SAW  DEVICE 


I  N1.  Alsup 

Naval  Undersea  C  enter 
San  Diego.  (  A  921 32 


ABSTRACT  A  metlunl  of  calculating  the  discrete  Fourier  transform  through  the  prime  transform  algorithm  with  surface 
acoustic  -.vase  device*  is  presented.  The  method  is  similar  to  the  chirp-Z  transform  (CZT)  technique  |4| .  and  utilizes  the  SAW  de¬ 
vice  as  a  transversal  I  liter  The  prune  transform  is  bawd  on  an  algorithm  which  uses  as  indices  the  set  of  integers  generated  by  the 
series  R  H  nodulo  N  (•  *  1 . 2. .  .  N  -  I )  N  is  prime  and  R  is  an  integer  whose  special  property  relative  to  N  is  that  its  successive 

powers  modulo  N  are  distinct,  and  thereform  form  a  permutation  of  the  integers  1.2 . N  -  I .  The  SAW  prime  transform  im¬ 

plementation  ha'  the  same  processing  speed  advantage  as  the  SAW  CZT  implementation,  namely,  that  it  computes  a  discrete  Fourier 
transform  with  speed  commensurate  to  a  fully  pipelined  FFT  running  at  the  same  sample  rate.  The  attributes  of  small  size,  light 
weight,  and  interconnection  simplicity  are  also  maintained. 


1  lie  Prime  1'ransform 

The  discrete  Fourier  transform  i  DFT)  of  a  sampled  data  se 
quencc  has  special  properties  when  the  number  of  points  to  be  trans 
formed  is  prime  1 1 1  For  each  prime  number  N  there  exist  integers, 
known  as  “primitive  roots.”  whose  successive  integer  powers  modulo 
N  will  generate  a  permuted  version  of  the  sequence  1 .  2. ....  N  -  I  |2J. 
The  DFT  <  I )  can  then  be  written  in  terns  of  the  permuted  integer  se¬ 
quence  for  neiwero  values  of  the  time  and  frequency  indices,  so  that 
(2)  results. 

N-l 

gncxp(-j:»mn/N>  (I) 
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»R.«P H>Kh+,/N>  «) 

C*l 

where 


m  =  Rk  mod  N.  k  r  0.  1 . 

, . .  N-2. 

n  *  Rk  mod  Y  f  •  1. 2. 

.N-l. 

Direct  permutation  of  sampled  Jiialog  data  can  be  achieved, 
for  example,  with  a  modified  version  of  a  commercially  jvai’  de¬ 
vice.'  flic  auxiliary  computation  to  account  lor  the  contribution  of 
the  zeroth  sample  can  be  carried  out  in  a  number  of  different  ways, 
as  are  shown  in  Figure  2. 


and 


f  igure  2  Four  Ways  to  Include  Auxiliary  Computation  for  Zeroth- 
Sample  Input 


Equation  (2)  represents  a  circular  convolution  with  auxiliary 
operations,  and  as  such  is  suitable  lor  implementation  by  tiansversal 
filter  architectures  (3)  The  filter  tap  weights  are  given  by  tile  expres 
sion  exp  (-j2»Rk/N).  and  are  permuted  values  of  a  complex  sinusoid 
The  sequence  g^j  is  a  permuted  version  of  the  input  sample  values 
with  the  first  sample  deleted  The  transform  output  coefficients 
are  the  m  *  I  to  N  -  1  Fourier  coefficients  in  permuted  order  and  Gc 
is  computed  separately  A  comparison  of  the  basic  prime  transform 
implementation  with  that  of  the  CZT  is  shown  in  Figure  I  (also  see 
14.51) 


»  K*iifv  I  mock  inaerum  t  ompurixon  i simplified!  between  Prime 
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For  a  given  pnme  N.  there  exists  a  selection  of  integers  between 
0  and  N  which  qualify  as  primitive  roots  The  number  of  such  roots  is 
given  by  the  Euler  ^-function.  OtN-1 ).  For  example,  there  arc  8  primi¬ 
tive  roots  associated  with  N  *  31  3.  1 1 .  I  2.  1 3.  1 7,  21 .  22.  and  24 

Thus,  there  are  (RN-I )  different  ways  to  fabricate  a  pnme  transform 
device  of  length  N  The  variety  of  related  pulse  compression  codes  is 
discussed  in  |7) .  Primitive  roots  for  various  primes  are  listed  in  181 
Equation  (3)  illustrates  the  Fuler  0-function 


0(A) 


a  -  I  .  b  -  1  c  -  1 
a  b  c 


»3) 


where  A  =  a°b^c7' 
are  integers 


and  a.  b.  c.  .  are  distinct  primes,  and  a.  *3.  ■>. 


1  Hr i icon  Corp  Mod-1  SAM  («mpled  analog  memory). 


•To  ht  pubtltM  In  Ultimonict  Sympottum  hoctr*ngt.  IEEE  (1975) 
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S  \W  Implementation 


\  surface  ucomti.  wave  (SAW)  device  was  designed  and  con¬ 
structed  j\  MX  to  ventv  the  feasibilit>  of  the  prime  transform  so 
implemented.  The  prime  number  was  chosen  to  be  31 .  and  the  primi¬ 
tive  root  3  This  led  to  the  selection  of  the  tap  weights  tabulated  in 
Table  I  Since  circular  convolution  was  required,  two  periods  less  one 
sample  ot  the  30-tap  sequence  were  incorporated  into  the  tapping 
structure,  resulting  in  a  total  of  59  complex  taps  The  complex  arith¬ 
metic  was  implemented  by  using  real  and  imaginary  parts  of  the  sped* 
fied  tap  response  to  determine  amplitude  weightings  on  two  parallel 
acoustic  paths  driven  b>  a  common  acoustic  input  signal. 


fable  I  Prime  Transform  Tap  Weights,  N  -  3 1 .  K  =  3 
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The  entire  electrode  structure  and  a  closup  are  shown  in  Fig¬ 
ures  3  and  4.  The  device  was  built  using  AT*cut  quart/  as  the  sub¬ 
strate.  with  tap  spacines  of  0256  inch,  equivalent  to  0.2  psec  Each 
tap  consists  of  8  finger-pairs  with  wavelength  0016  inch,  equivalent 
to  80  MHz.  The  taps  are  spaced  with  a  gap  equal  to  their  width,  and 
the  entire  array  of  taps  with  one  input  transducer  is  1 .88  inches  in 
length  flic  aperture  of  the  taps  varies  according  to  the  specified 
weighting  function,  and  has  a  maximum  of  0  256  inch  The  aperture 
ot  the  input  transducer  at  either  end  of  the  device  is  0  596  in.h  The 
interconnect  busses  were  configured  so  that  either  all  59  or  just  30 
taps  (from  either  end)  could  be  connected  together  to  form  the  real 
and  imaginary  output  signal  components 


I  igtirc  3  Prune  Ir.mstmm  S \W  Electrode  Structure 


Figure  4.  Prime  Transform  SAW  Electrode  Structure.  Closeup 


A  separate  uniformly  weighted  tapping  structure  (not  shown) 

I  was  designed  to  carry  out  the  /eroth-sample  auxiliary  calculation,  but 
was  not  used  in  initial  device  tests  The  resulting  device  is  capable  of 
calculating  a  31  point  DFT  in  6.2  psec.  and  can  operate  at  duty  cycles 
up  to  50 'y  Two  such  devices  can  be  used  alternately  to  achieve  100'* 
duty  cycle  when  required 

Experimental  Results 

The  procedure  used  to  view  the  impulse  response  of  the  S  \W 
prime  transform  device  is  illustrated  in  Figure  5  A  quarter-wave  delay 
(about  3.1  ns)  was  used  to  enable  quadrature  combination  of  the  real 
and  imaginary  parts  of  the  response,  and  the  separate  and  combined 
signals  arc  shown  in  Figures  6  and  "*  I  his  technique  for  achieving  dis¬ 
crete  phase  modulation  of  a  signal  earner  is  discussed  in  |6| .  and 
could  have  just  as  easily  been  incorporated  into  the  initial  mask  de¬ 
sign  instead  of  being  left  to  electrical  manipulation  external  to  the 
SAW  device  as  performed  here  I  jck  of  perfection  in  implementa¬ 
tion  of  the  spec  died  tap  weights  due  to  missing  or  broken  fingers  will 
show  up  as  a  nonuniformity  in  the  complex  output  sample  magnitude, 
and  f  igure  7  does  show  this  to  some  extent 
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Figure  "a  Impulse  Recponse.  Tips  30-1.  Complex  Combination 


Figure  9a  S  \\V  Prime  Transform  Outputs  at  Carrier  <  SO  MHz)  Upper. 
Input  frequency  =  I  cycle.  Lower.  Input  Frequency  * 

30  cycles. 


Figure  10  SAW  Prime  Transform  Outputs,  Magnitude  Squared,  at 
Baseband:  Upper.  Input  Frequency  *  I  cycle;  Lower, 
Input  Frequency  =  30  cycles. 

Conclusions 

The  prime  transform  algorithm  implemented  with  SAW  trans¬ 
versal  filters  has  been  shown  to  be  a  powerful  tool  for  the  high-speed 
evaluation  of  the  discrete  Fourier  transform  At  a  5  MHz  sample  rate, 
a  3 1 -point  OFT  can  be  calculated  in  6.2  //sec.  which  represents  a  sub¬ 
stantial  increase  in  speed  over  conventional  FIT  techniques.  The 
transversal  filter  prime  transform  is  very  similar  to  the  CZT  method, 
and  represents  an  alternative  to  it  when  permutation  of  the  data  being 
transformed  is  preferred  over  multiplication  of  the  data.  The  prime 
transform  requires  auxiliary  calculations  to  account  for  the  contribu¬ 
tion  of  the  zeroth  data  sample  and  to  evaluate  the  dc  transform  coef¬ 
ficient.  but  either  or  both  of  these  may  be  dispensed  with  if  the  ap¬ 
plication  does  not  depend  critically  upon  them  Inverse  permutation 
of  the  filter  output  may  also  be  eliminated  if  not  needed. 

The  same  limitations  with  regard  to  accuracy  apply  to  analog 
implementations  of  the  prime  transform  algorithm  such  as  this  as  also 
apply  to  CZT  analog  implementations,  which  limitation  is  about  the 
equivalent  of  7  or  8  bits.  Like  the  CZT  filter,  the  transversal  filter 
prime  transform  implementation  affords  the  possibility  of  calculating 
both  forward  and  inverse  transforms  using  identical  modules,  and 
this  in  turn  leads  to  the  capability  of  real-time  high-data-rate  linear 
signal  processing  tools  such  as  cross-convolvers,  cross-correlators, 
adaptive  filters,  and  programmable  matched  filters. 
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J.M.  Speiser 
Naval  Undersea  Center 

ABSTRACT 

A  modular  architecture  is  described  for  the  implementation  of  the  Even 
Discrete  Cosine  Transform  (EDCT).  This  architecture  permits  the  use  of  four  EDCT 
modules  to  compute  a  double  length  transform  with  twice  the  throughput  rate  of 
the  individual  modules. 

INTRODUCTION 

The  utility  of  the  EDCT  for  data  compression  has  been  previously 
described  [1-3].  Several  different  serial  access  implementations  for  its  high 
speed  computation  have  also  been  described  [4].  Implementations  using  charge 
coupled  device  (CCD)  transversal  filters  are  particularly  attractive  [2]  for 
applications  which  require  low  weight,  small  size,  low  power  consumption,  and 
controllable  clocking  of  the  cpmputation.  Present  CCD  transversal  filters  perform 
well  at  shift  rates  of  up  to  about  5  x  10^  samples  per  second.  This  is  a  factor 
of  two  too  slow  to  handle  conventional  television  signals  using  the  previously 
described  EDCT  architectures.  Thisnote  describes  a  subdivision  of  the  computation 
tasks  to  permit  greater  parallelism  in  the  hardware  to  increase  both  the  throughput 
and  transform  size  implementable  with  a  fixed  set  of  transversal  filter  and  chirp 
read-only  memory  modules. 

DCT  IMPLEMENTATIONS 

A  system  to  compute  the  EDCT  of  length  N,  implementing  previously 
given  equations  [4]  is  shown  in  Fig.  1.  It  differs  from  serial  access  implementations 
of  the  ODCT  [2]  primarily  in  the  sinusoidal  multiplication  following  the  chirp 
postmultiplication.  Despite  this  slight  complication  the  EDCT  was  selected  for 
modular  decomposition  rather  than  the  ODCT  because  an  odd  length  extension  of  a 
data  block  can  only  be  subdivided  into  a  even  length  data  block  and  an  odd  length 
data  block,  while  an  even  length  extension  of  a  data  block  can  be  subdivided  into 
two  data  bit  xs  of  the  same  length  in  order  to  permit  simultaneous  computation 
by  similar  modules. 

Since  intermediate  complex  quantities  need  to  be  preserved  in  the  modular 
decomposition,  the  complex  extended  discrete  Fourier  transform  (EDFT)  portion  of 


Fig.  1  was  chosen  as  the  basic  module.  This  module  computes  N  points  of  a  length 
2N  DFT  of  the  extension  of  the  data  sequence  by  N  zeroes.  The  interconnection  of 
modules  with  minor  auxilliary  components  to  perform  an  EDCT  with  doubled  throughput 
on  a  double  length  data  block  is  shown  in  Fig.  2. 


DERIVATION  OF  THE  MODULAR  EDCT 

Let  the  input  data  be  denoted  by  gm  for  m=  0,1, ...M-l.  Define  the 
symmetrized  extension  of  the  data  by  g_1_m  =  gm  for  m=0,l , . . .M-l .  The  even 
discrete  cosine  transform  (EDCT)  of  g  may  then  be  defined  by  any  of  the  equivalent 
expressions  shown  in  equations  1-3. 
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Let  the  data  block  length  be  even,  say  M=2N.  Then  the  summation  in 
the  DCT  can  be  split  into  shorter  sums  as  shown  in  equations  4-5. 
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Each  of  the  summations  in  equation  5  may  be  interpreted  as  a  OFT  of  length 
2N  of  a  length  N  data  block  extended  by  N  zeroes.  The  first  N  coefficients  of  such 
a  OFT  are  computed  directly  by  the  EDFT  module  of  Fig.  1.  The  remaining  coefficients 
may  all  be  written  in  :the  form  G^,  where  p=0,  1....N-1.  The  corresponding  terms 
of  equation  5  are  examined  in  equations  6-9. 
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The  structure  shown  in  Figure  2  uses  equation  4  to  generate  the 
coefficients  Gq,G^ , . . .G^_^  and  uses  equations  6-9  to  generate  GN»G^+^, - ^N-l 
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MODULE  TO  PERFORM  OFT  OF  SIZE  2N  OF  LENGTH  N  DATA 
EXTENDED  BY  N  ZEROS  (EDFT) 


SYSTEM  TO  PERFORM  AN  EDCT  OF  LENGTH  N 


System  to  Perforin  an  F.DCT  of  Length  N  Including  Module  to  Perforin  a  DFT  of  Size  IN 
of  a  Data  Block  of  Size  N  Fx  tended  by  N  Zeros. 


89 


#*«*>*» 


APPENDIX  G 

SIGNAL  PROCESSING  ARCHITECTURES 
USING  TRANSVERSAL  FILTER  TECHNOLOGY 


I 


SIGNAL  PROCESSING  ARCHITECTURES  USING 
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INTRODUCTION 

A  large  portion  of  the  computational  load  for  many  signal  processing  problems  consists  of  the  com¬ 
putation  of  linear  transforms  Lor  time-invariant  linear  transforms  such  as  cross  convolution  or  matched 
filtering,  the  transversal  filter  provides  a  highly  parallel  computational  module  with  high  throughput  and 
minimal  control  overhead  ( I  | .  This  paper  will  show  how  similar  computational  modules  can  he  configured 
to  provide  similar  computational  advantages  for  a  large  class  of  time-variant  linear  transforms  including 
one-dimensional  and  multi-dimensional  discrete  Fourier  transforms  and  one-dimensional  and  two- 
dimensional  discrete  cosine  transforms.  Furthermore,  time-variant  transform  modules  may  he  combined 
to  implement  high  capacity  time-invariant  linear  transforms.  The  implementation  of  these  techniques  using 
surface  acoustic  wave  (SAW)  and  charge  coupled  device  (CCD)  technology  permits  the  real-time  solution  of 
several  important  signal  processing  problems  including  image  data  compression,  wideband  radar  signal 
analysis  and  spread  spectrum  communications. 


I 


l 

I 


COMPUTATIONAL  MODULES 

A  linear  transform  on  sampled  data  of  finite  extent  may  be  viewed  as  the  multiplication  of  a  vec¬ 
tor  by  a  matrix.  Multiplication  by  diagonal,  circulant.  or  Toeplitz  matrices  may  be  accomplished  rapidly 
with  simple  computational  hardware  modules.  Multiplication  by  an  N  X  N  diagonal  matrix  requires  only 
a  scalar  multiplier  and  a  memory  containing  N  values  to  provide  serial  access  to  the  reference  function 
Multiplication  by  an  N  X  N  Toeplitz  matrix  corresponds  to  a  convolution  and  may  be  performed  using 
a  transversal  filter  having  2N-I  taps  Multiplication  by  an  N  X  N  circulant  is  a  special  case  of  multiplica¬ 
tion  by  N  X  N  Toeplitz  matrix  in  which  the  length  of  the  transversal  filter  may  be  reduced  to  N  taps  it 
the  data  block  is  recirculated  through  the  filter  or  reread  into  the  filter  from  a  buffer  memory. 


ONE-DIMENSIONAL  DFT 

Linear  filters  have  been  used  for  many  years  for  the  calculation  of  the  power  spectra  of  continuous 
signals.  One  of  the  earliest  methods  used  a  bank  of  wave  filters  to  measure  the  spectra  in  fractional  octave 
bands  for  telephone  network  equalization  (  21 .  However,  when  increased  resolution  was  required  the  num¬ 
ber  of  filters  rapidly  become  unmanageable.  An  alternative  which  overcame  the  difficulty  of  a  large  num¬ 
ber  of  filters  each  with  small  time-bandwidth  product  was  to  substitute  one  linear  I'm  (chirp)  filter  with 
large  time-bandwidth  product  and  to  employ  matched  filtering.  In  this  system  the  signal  to  be  analyzed 
is  used  to  single  sideband  (SSB)  modulate  a  locally  generated  chirp  signal  and  the  composite  modulated 
signal  is  filtered  in  a  chirp  delay  line  matched  filter.  Each  component  of  the  input  signal  spectrum  shifts 
the  locally  generated  chirp  to  a  different  position  in  the  spectrum  after  SSB  modulation  and  these  shifted 
chirps  then  correlate  with  the  reference  signal  represented  as  the  impulse  response  of  the  matched  filter  at 
different  times.  Thus  the  output  signal  amplitude-time  history  reflects  the  amplitude-frequency  composi¬ 
tion  of  the  input  signal. 


•To  be  published  in  Proceedings  of  the  19  75  IEEE  International  Symposium  on  Circuits  and  Systems  (Boston.  21-23  April.  1975). 


93 


Bleustein  |3|  recognized  (hat  the  discrete  homier  transform  (Oh  I  )  ol  sampled  data  was  amenable 
to  a  similar  interpretation.  In  addition  to  just  calculating  the  magnitude  of  the  Fourier  transform,  linear 
filters  could  calculate  the  phase  and  thus  all  ot  the  operations  such  a  cross  convolution  and  a  cross  correla¬ 
tion  could  he  calculated.  This  technique  came  to  be  called  the  chirp-/,  transform  <(  Zl  )  and  can  he  applied 
to  other  problems  besides  just  the  calculation  of  the  DFT  |4| .  Prior  fo  these  developments,  digital  compu¬ 
tation  of  the  DFT  had  been  significantly  improved  by  the  use  of  a  special  algorithm  called  flic  last  Fourier 
transform  ( I- FT )  which  was  described  by  Cooley  and  Tukey  |5  I .  The  FFT  algorithm  gained  rapid  popu¬ 
larity  in  signal  processing  since  it  allowed  the  calculation  ol  the  DF  T  to  be  done  using  significantly  fewer 
machine  operations  ( multiplications)  than  direct  evaluation. 

The  DFT  can  be  defined  with  many  normalizations  just  as  the  continuous  Fourier  transform  can 
be.  For  this  paper  the  definition  of  the  direct  or  forward  transform  of  a  complex  vector  g  of  length  N 
is  given  as 


N-l 

G'  =  F*  =  £  e-J’r-nm/N  gn 
n=0 

and  the  inverse  transform  as 

n=(J 


where  F  is  an  N  X  N  matrix  with  elements  P'n  m  =  ^ 

By  direct  inspection  it  is  observed  that,  if  symmetries  of  the  function  expjffJnm/N  are  not  ex¬ 
ploited.  then  the  number  of  complex  multiplications  required  will  be  N~  corresponding  to  N  multiplica¬ 
tions  for  each  frequency  component  evaluated.  Fven  on  high  speed  digital  computers  this  can  become 
the  limiting  consideration  in  signal  processing  applications.  The  advantage  of  the  FFT  algorithm  is  that 
for  highly  composite  values  of  the  DFT  size  N  the  number  of  multiplications  is  proportional  to  N  logj  N. 

Although  the  FFT  has  been  successful  in  substantially  reducing  the  computing  time  and  cost  of 
using  general  purpose  digital  computers  it  has  several  disadvantages  lor  special  purpose  real  time  computa¬ 
tion.  At  high  throughput  rates  which  are  required  for  real  time  image  processing  the  processor  either  must 
operate  logsN  times  faster  than  the  data  rate  or  pipeline  structures  which  use  distributed  memory  and 
log-iN  multipliers  must  be  used.  In  addition,  the  internal  arithmetic  ol  the  III  processor  must  be  done  at 
increased  precision  in  order  to  compensate  for  the  multiple  round  oil  errors  introduced  by  the  successive 
stages  in  the  FFT  processor.  Although  these  difficulties  can  be  overcome,  it  is  not  always  possible  to 
arrange  the  computation  in  a  form  where  the  size  of  the  transform  is  highly  composite.  For  the  above  rea¬ 
sons  and  because  of  the  difficulty  of  obtaining  small,  low  power,  last  analog  to  digital  converters,  linear 
transversal  filter  inplementations  of  the  chirp-Z-transform  are  attractive  1 0 1  rather  than  the  previous  (  ZT 
implementation  which  used  an  FFT  to  perform  the  required  convolution. 

The  DFT  may  be  easily  reduced  to  the  form  suitable  for  linear  filtering  by  the  substitution 


which  changes  a  product  of  variables  into  a  difference  so  that 


N- 


r;  =  ,.-jtrm-  N  V  Jjr<n-m>_  N  -j7rn-  N  „ 
'’m  c  ^  *■  1  hi 

n=U 


<  4 1 


94 


I  Iiin  lorin  is  seen  In  be  k*L|iti v.iK-n l  in  factoring  the  Fourier  matrix  F  into  the  product  ot  three  matrue' 


F  DID  i'. 

■) 

where  I)  o  a  di.ieon.il  m.itiix  with  elements d|jn  =  exp  ( —  j tt n _  N  l  and  1  is  a  I oeplit/  matrix  with  elements 
l'nin  =  exp  (tirtn  -ml  Ni 

I  he  t  /  I  .ilgoiithtn  is  easily  implemented  by  transversal  I  liter  techniques  In  this  ease  the  1)1  I 
is  computed  by  piemultiplnation  by  a  discrete  chirp,  convolution  with  a  discrete  chirp,  and  postmultipli- 
cation  In  a  discrete  cluip  I  igure  I  shows  this  configuration.  However,  it  must  he  remembered  that  both 
the  multiplications  and  convolutions  are  complex  and  a  suitable  representation  of  the  complex  numbers 
must  be  used.  One  representation  is  by  real  and  imaginary  part.  Figure  2  shows  the  DFT  organized  as  a 
CZT  and  implemented  with  parallel  computation  of  the  real  and  imaginary  parts.  In  Figure  2  the  input 
signal  is  represented  as  g  =  g^  +  jgj  and  the  output  signal  is  represented  as  (,'  =  +  j(/.  where  it  is  under¬ 
stood  that  g  =  -a,,  n  =  0.  .  N  -  I  and  G  =  C»n  n  =  0 . N  -  I 

In  order  to  determine  the  specific  form  of  the  transversal  filters  it  is  necessary  to  know  the  specific 
value  of  N.  When  N  is  odd  the  Toeplitz  matrix  T  may  be  represented  as  a  transversal  filter  with  2N  -  I  com¬ 
plex  tapsh_(N_j  ^  to  hjyj_j  where  hn  =  W_n  n  =  -<N-h  to  (N-l  I.  and  W  =  exp  <-j2ir/N>.  The  required 
convolution  has  been  implemented  with  the  general  transversal  filter  shown  in  Figure  3. 

When  N  is  even,  it  can  be  shown  that  Tn  m  =  T^+n  m  where  the  subscripts  are  reduced  mod  N. 

rims  I  is  a  circulant  matrix  and  can  be  implemented  with  a  recirculating  transversal  filter  as  shown  in 

*> 

Figure  4  where  the  number  of  complex  taps  is  N  and  tap  weights  are:  hn  -  n  =  0 . N  -  I . 

In  tile  specific  case  when  N  is  an  odd  prime,  additional  simplification  is  possible.  It  is  possible  to 
eliminate  the  multipliers.  The  DFT  may  be  written  as 


N=1 


^o  ~ 

£  h 

(0) 

n=() 

^’m  " 

h  -  E 

for  m  =  1 .  . 

.  .  .  N  -  1 

(71 

n=  I 


*  Denotes  either  convolution  or  circular  convolution 


Figure  I  (  hirp-Z-Transform  Implementation  of  the  DFT 
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f  igure  4.  Circular  Convolution 


Since  only  nonzero  values  of  n  and  m  occur  in  the  right  hand  side  of  fcquation  (7)  it  is  possible  to 
use  Gauss’s  analogy  between  logarithms  and  indices  with  respect  to  a  primitive  root  |  7 1  to  replace  the 
product  nm  by  a  primitive  root  raised  to  a  sum.  thus  reducing  this  computation  to  a  correlation  between 
permuted  functions  1 8 1 .  In  matrix  notation 

(!'  -  gQ  !  =  F'  g'  (8a) 


where 


/  = 


(Xb) 


and  the  matrix  F'  can  be  factored  into  the  three  matrices 


F'  P*CP. 


where  P  is  a  N  -  I  permutation  matrix.  C  is  a  ( N  -  I  X  N  -  I )  circulant  matrix.  P1  is  the  transpose  of  P.  (!' 
and  g'  are  column  vectors  of  size  N  -  I  derived  from  O'  and  g  by  deleting,  respectively.  CQ  and  gQ  The 
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elements  of  the  matrix  (' are  C|i  ni  =  F(  n)  (m)  n.  in  -  I . N  -  I  and  <n)p  -  r11  mod  N  .  w  here  i  i^  a 

primitive  root  of  N  and  the  elements  ot  P  are 


I)  =  x 
1  nm  "(nlptn 


n.  m  =  I _ N  -  I 


(10) 


where  6^,  is  the  Kronecher  delta  and  p  implies  permuted. 

Thus,  for  the  ease  where  N  is  an  odd  prime.  C  is  a  circulant  matrix  and  can  he  implemented  with  a 
recirculating  transversal  filter  whose  tap  weights  are 


II 

P 

n  =  1  . 

.  .  N  -  1 

(  1  la) 

(n)n 
=  W  P 

n  =  1 . 

.  .  .  N  -  1 

(lib) 

and  the  Fourier  transform  coefficients  are  given  by 


N-l 

G0  ~  ^  1  Sn 
n=0 


( 1  2a  > 


and 


G’m  Gm  +  gQ 


m  =  I  .  .  .  N  -  1 


(12b) 


These  concepts  are  illustrated  in  Figure  5  for  N  =  5.  Thus,  the  one-dimensional  architectures  may  be  sum¬ 
marized  as  a  time-varying  operation,  a  convolution,  and  a  second  time-varying  operation  where  for  arbitrary 
N  the  time-varying  operations  are  multiplication  by  a  diagonal  matrix:  and  where  for  the  special  case.  N  is 
an  odd  prime,  the  time-varying  operations  are  multiplication  by  a  permutation  matrix  and  its  transpose. 


Figure  5.  I  xample  of  the  Prune  Transform  when  N  =  5 
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IMPLEMENTATION 


Many  types  of  transversal  lillei  implementations  may  he  easily  accomplished  using  analog  sampled 
data  techniques  |  v| .  Each  oilers  some  advantages  and  some  disadvantages  in  any  particular  application, 
although  they  are  all  architecturally  similar.  In  this  paper  the  implementations  to  he  discussed  are  surlace 
acoustic  wave  (SAW)  devices,  charge  coupled  devices  (CCD)  and  hybrid  analog-digital  correlators. 


SURFACE  ACOUSTIC  WAVE  DEVICES 

Surlace  Acoustic  Wave  (SAW)  devices  can  accept  either  analog  or  sampled  analog  input  and  the  out¬ 
put  is  analog.  A  substrate  of  piezoelectric  material  is  polished  on  one  face  and  a  pattern  ol  aluminum  or 
other  conductor  is  deposited  by  photolithographic  techniques.  In  its  simplest  configuration,  sets  ot  inter 
digitated  finger  electrodes  are  spaced  at  the  sampling  rate  distance  through  the  use  of  the  relationship  d  - 
Co  t  where  d,  is  the  tap  spacing.  CK  is  the  Rayleigh  wave  velocity,  and  ts  is  the  sampling  increment  l  or 
a  typical  substrate  of  ST-eut  quart/. CR  ~  3mm/psee  and  for  a  typical  sampling  increment  ol  ISO  nsec, 
d  =  0  45  mm  and  a  1 50  point  Dl'f  can  be  implemented  by  recirculating  convolution  in  an  active  length  less 
than  75  mm  (3  in.)  on  a  100  mm  (4  in.)  substrate.  A  simple  real  transversal  filter  is  shown  in  Figure  h  along 
with  the  method  of  establishing  the  lap  weights.  A  prototype  complex  filter  1 101  is  shown  packaged  in 
Figure  7.  A  similar  device  has  been  used  as  the  premultiplication  and  postmultiplication  reference  function 
generator  and  double-balanced  mixers  have  been  employed  as  the  multipliers 

Physical  limitations  on  the  si/e  of  available  substrates  limit  the  available  DFT  size.  For  an  active 
length  of  substrate  L  =  150  nun  and  a  surface  wave  velocity  CR  =  3  mm  '/zsec  the  maximum  transform 
size  Nmax  for  a  data  rate  Fs  is  approximately  Nmax  =  LSFS/CR  -  50  Fs  with  Fs  in  megaHertz  Current  SAW 

111  d  A  •’ 


Figure  6.  Surface  Wave  Transversal  Filter 
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C2T  TRANSVERSAL  FILTER 


Figure  7.  32  Complex  Tap  SAW  CZT  Filter 


transversal  filters  operate  at  sample  rates  from  1  to  100  MHz  with  carrier  frequencies  from  5  to  500  MFlz 
with  typical  fractional  bandwidths  of  about  I  O'-' .  Although  there  are  some  piezoelectric  semiconductors, 
silicon  is  nonpiezoelectric.  Heteroepitaxial  materials  such  as  Aluminum  Nitride  (AIN  )  and  Silicon  (Si) 
cogrown  on  a  sapphire  substrate  are  currently  under  development  for  monolithic  integrated  circuits  where 
The  AIN  is  the  piezoelectric  and  Si  is  the  semiconductor.  Such  monolithic  circuits  should  make  possible 
monolithic  CZTs  at  100  to  500  MHz  data  rates. 


CHARGE  COUPLED  DEVICES 

CCDs  are  sampled  data  analog  circuits  which  can  be  fabricated  by  Metal  Oxide  Semiconductor 
(MOS)  technology  as  LSI  components  |  I  1 1 .  As  such  they  are  directly  compatible  with  other  MOS  circuits 
Current  CCD  transversal  filters  have  operated  as  video  devices  with  sample  rates  up  to  5  MHz.  CCDs  oper¬ 
ate  by  the  manipulation  of  injected  minority  carriers  in  potential  wells  under  MOS  capacitors  and  thus 
behave  as  capacitive  reactances  with  low  power  dissipation.  However,  since  the  potential  wells  which  con¬ 
tain  the  minority  carriers  also  attract  thermally  generated  minority  carriers,  there  is  a  maximum  storage 
time  for  the  analog  signal  which  depends  on  the  dark  current  associated  with  the  temperature  of  the  silicon 
Under  normal  conditions  at  room  temperature,  dark  currents  are  tens  of  nAmps/em“  and  storage  times  of 
hundreds  of  milliseconds  can  be  achieved. 

There  are  many  ways  in  which  unidirectional  charge  transfer  can  be  achieved  The  first  developed 
was  a  three-phase  clocking  structure  which  is  illustrated  in  the  transversal  filter  of  Figure  b.  The  three 
electrode  CCD  structure  is  planar,  much  like  the  SAW  devices,  and  the  direction  of  charge  propagation  is 
determined  by  the  sequence  of  potentials  applied  to  the  three  electrodes.  Unfortunately,  if  the  minority 
carrier'  are  allowed  to  collect  at  the  semiconductor-oxide  boundary,  poor  charge  transfer  efficiency  will 
result  due  to  minority  carriers  getting  caught  in  trapping  sites.  Fins  means  that  the  CCD  will  behave 
nonlinearly  unless  there  is  sufficient  propagating  charge  present  to  fill  all  of  the  traps  By  biasing  the 
operating  condition  of  the  CCD  so  that  about  IOC  of  the  dynamic  range  is  used  for  the  injection  of  a  "tat 
zero."  the  traps  are  kept  continuou  ly  tilled  and  the  device  It. is  over  a  Ml  dB  dynamic  range  In  practice. 


I • » t*u rc  N  Schematic  of  the  Sampling  W’oijih I niji.  atnl  Summing  Operation 


a  video  signal  representing  the  signal  to  he  processed  is  added  to  a  fixed  bias  somewhat  larger  than  one-halt 
of  the  peak-to-peak  value  of  the  signal  Since  the  effective  storage  tune  of  the  device  is  long  relative  to  the 
time  required  to  execute  a  convolution.  CCDs  can  he  considered  to  he  interruptible  signal  processors  and 
as  such  are  more  compatible  with  the  executive  control  required  for  signal  processing.  A  h4  point  CCD 
lilter  with  discrete  cosine  transform  sine  and  cosine  chirps  is  shown  in  Figure  0.  This  chip  was  developed 
bv  1  exas  Instruments  for  the  Naval  Undersea  Center  for  image  processing  The  discrete  cosine  transform  is 
described  in  a  subsequent  section. 

Current  research  in  CCDs  is  directed  toward  improving  the  charge  transfer  etticiency  and  removing 
the  requirement  of  continuous  “fat  zero"  charge  injection  by  ion  implantation  techniques  which  keep  the 
minority  carriers  away  from  the  semiconductor  oxide  boundary  Ion  implantation  is  also  being  used  to 
provide  asymmetric  potential  wells  so  that  simpler  two-phase  clocking  can  he  employed.  Currently  avail¬ 
able  CCDs  have  500  stages  with  0  0000  transfer  efficiency  and  devices  w  ith  up  to  2000  stages  are  planned 

Another  charge  transfer  device  similar  to  the  CCD  is  the  Bucket  Brigade  Device  (BBI)i  I  his  is  a 
sequence  of  MOS  transistors  coupled  together  by  diffusion  enhanced  Miller  capicitance.  Although  these 
devices  do  not  operate  at  frequencies  as  high  as  CCDs,  they  have  better  low  trequency  pertornianee  since 
t hey  include  active  devices.  A  CZT  has  been  implemented  with  two  BBI)  chips.  I  wo  200  tap  filters  are 
implemented  on  each  chip:  one  a  discrete  cosine  and  the  other  a  discrete  sine  lilter  1  he  device,  the 
complex  chirp  used  in  the  premultiplier  and  a  typical  input  and  output  are  shown  in  F  igure  10.  The  in¬ 
put  is  an  offset  cosine  wave  and  the  output  shows  a  D  C  component  plus  a  response  at  the  cosine  wave 
trequency.  These  filters  can  operate  at  100  k  1 1  /  and  have  tap  accuracies  bcttei  than  |\  Skith  careful  con¬ 
trol  of  geometry,  both  BBD  and  CCD  filters  with  tap  accuracies  approaching  0  I  should  be  possible  This 
chip  was  also  developed  by  Texas  Instruments  for  the  Naval  Undersea  (  enter 
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Figure  l).  o4  Point  CCD  Filters 


A  hybrid  binary-analog  CCD  correlator  module  with  32  analog  taps  1 12|  has  been  made  for  the 
Naval  Undersea  Center  by  General  Hlectric  and  is  shown  diagrammatically  in  Figure  I  I  The  maximum 
clock  rate  of  the  module  is  4  MHz.  These  modules  may  be  cascaded  to  increase  the  correlation  length  to 
provide  an  analog  versus  multilevel  cross  correlator  with  analog  output.  The  correlator  module  uses  charge 
propagation  through  only  3  stages  since  the  CCD  is  arranged  in  32  stages  of  3  samples  each.  By  transverse 
shifting  of  the  charge  instead  of  longitudinal  shifting,  charge  transfer  inefficiency  degradation  is  avoided 
and  the  modules  can  be  cascaded  to  very  large  sizes.  However,  with  this  configuration,  dark  current  non¬ 
uniformity  must  be  controlled  to  prevent  time  variable  pattern  noise  from  contributing  to  the  output.  Sys¬ 
tematic  errors  in  the  tap  weights  can  be  measured  and  stored  and  the  analog  signal  corrected  for  the  mea¬ 
sured  variation  in  uniformity  before  it  is  stored  in  the  CCD  registers. 


DISCRETE  COSINE  TRANSFORM 

Closely  related  to  the  DFT  is  the  discrete  cosine  transform  (IX'T).  Two  different  types  of  DCTs 
are  useful  for  reduced  redundancy  television  image  transmission.  Both  are  obtained  by  extending  a  length 
\  real  data  Nock  to  have  even  symmetry,  taking  the  discrete  Fourier  transform  1 1)1- T)  of  the  extended  data 
block,  and  saving  N  terms  of  the  resulting  DFT  Since  the  DFT  of  a  real,  even  sequence,  is  a  real,  even  se¬ 
quence.  either  IX'T  is  its  own  inverse  if  a  normali/cd  1)1  I  is  used. 

The  "Odd  DCT"  (ODCT)  extends  the  length  N  data  block  to  length  2N-I .  with  the  middle  point 
of  the  extended  block  as  a  center  of  even  symmetry.  The  "Even  DCF"  (I  DCT)  extends  the  length  N  data 
block  to  length  2N.  with  a  center  of  even  symmetry  located  between  the  two  points  nearest  the  middle. 
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For  example,  the  odd  length  extension  of  the  sequence  g^g^g-,  is  gsg  |gQg|g->.  and  the  even  length  exten¬ 
sion  is  gig | gQgQg  ] g-> •  ln  both  cases,  the  symmetri/.ation  eliminates  the  jumps  in  the  periodic  extension  of 
the  data  block  which  would  occur  it  one  edge  of  the  data  block  had  a  high  value  and  the  other  edge  had  a 
low  value:  in  effect  it  performs  a  sort  of  smoothing  operation  with  no  loss  of  information.  It  will  be 
noted  that  the  terms  "odd”  and  “even"  in  OIKT  and  FDC'T  refer  only  to  the  length  of  the  extended  data 
block  in  both  cases  the  extended  data  block  has  even  symmetry.  Both  types  of  DCT  may  be  implemented 
using  compact,  high  speed,  serial  access  hardware,  in  structures  similar  to  those  previously  described  for 
the  chirp-Z  transform  (CZT)  implementation  of  the  DPT. 

Let  the  data  sequence  be  gQ.  gj ,  .  g^_j  The  OIX'T  of  g  is  defined  as 

N-l  -j-ffnk 

C,k  =  gn  e  2N-'  for  k  =  0.  I . N-l  (13) 

n=  -(N-l ) 


* 


|£<1-b  I  9 _ 


REFERENCE 

WORD 

BINARY 


,  REFERENCE 
‘shift  CLOCK 


Figure  I  1  Schematic  Diagram  of  a  Binary-Analog  CCD  Correlator  Module 

where  g_n  =  gn  for  n  =  0.  I ,  .  .  .  ,  N  -  I .  By  straight  forward  substitution  it  may  be  shown  that 
N-l  -j2trnk 

=  r„  e  -N_l 


where  g  is  defined  by  equation  ( 14a) 

_  _  j  0.5  gg,  n  =  0 

gn_  {  g„.n=l . N-l 


The  identity  (15)  may  be  used  to  obtain  the  CZT  form  of  the  ODCT  shown  in  equation  ( 16). 


2nk  =  n“  +  k“  -  (n  -  k)" 


-jffk-  N-l  -jffn- 


jtr(n-k)-  ) 


Gk  =  2Rej  e^"1  Z 


2N-I  ~ 
e  P.i  c 


Jtr(n-k 

2N-I 


I  he  I  DC  I  of  g  is  defined  by  equation  (  I  "'at.  where  the  extended  sequence  is  defined  In  equation 


i  I  ~  h ) 


(ik=e 


-jjrk_  n-I 

:n  £ 

n=-N 


-j27rnk 

2N 


for  k  =  0.  1 , 


g_l_n=gn  for  11  =  0.1 . N-l 

If  the  mutually  complex  conjugate  terms  in  equation  ( 17a)  are  combined,  then  equation  ( 18)  results 
Equation  (  18)  may  be  viewed  as  an  alternate  way  of  defining  the  LDC'T. 


■j2trnk 

2N 


I  -jrrk  n-I  : 
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Equation  ( 18)  may  be  put  in  the  C7.T  format  as 
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Ahmed  |  131  has  investigated  the  use  of  the  HIX'T  as  a  substitute  for  the  Karhunen-1  oeve  trans¬ 
form  for  exponentially  correlated  data  and  finds  that  it  is  superior  to  the  Fourier  transform  and  is  com¬ 
parable  to  the  Karhunen-Loeve  (K-L)  in  rate-distortion  performance  while  maintaining  the  computational 
simplicity  of  a  transform  which  does  not  depend  on  the  picture  statistics.  Habibi  1 14 1  has  shown  by  simu¬ 
lation  that  the  IX'T  is  equivalent  in  a  inean-square-error  sense  to  the  K-L  transform  under  basis  restriction. 
A  1 00  point  ODC'T  for  use  in  a  10  frame  per  second  experimental  TV  image  compression  system  has  been 
constructed  by  the  Naval  Undersea  Center  for  the  Advanced  Research  Projects  Agency  (ARPA)  using 
bucket  brigade  device  transversal  filters.  The  transversal  filter  has  199  nonzero  taps  and  both  the  cosine 
filter  and  sine  filter  required  for  complex  arithmetic  are  implemented  on  a  single  chip.  A  microphotograph 
of  the  filters  is  shown  in  Figure  1  2. 


TWO-DIMENSIONAL  DFT 

The  DFT  of  a  two-dimensional  array,  gp  ^  n  7  may  be  computed  by  successive  applications  of  the  one- 

dimensional  DFT.  This  concatenation  is  readily  observed  by  writing  the  expression  for  the  two-dimensional 
DFT  as 


Ns- 1  |  N , -I  ~i:7rnlkl  \  -j 27rn 2k 2 

(i  (k  | .  ki)  =  ^  v  ^  g  <  n  | .  n  2 )  e  N’  >  e  N- 

111=0  /  n  |  =0  l 
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Figure  12.  Bucket  Brigade  Chirp  Filters 


In  two  dimensions  the  concatenation  proper!  \  of  the  DFT  can  be  exploited  along  with  auxiliar>' 
memory  to  compute  the  two-dimensional  transform  by  successively  computing  the  CZT  of  rows  of  the 
input  signal  matrix  and  then  using  the  auxiliary  memory  as  a  row  to  column  transformation,  i.e.,  transpos¬ 
ing  the  partial  Fourier  transform  matrix,  and  computing  the  two-dimensional  Fourier  transform  with  a 
second  CZT. 

An  alternative  method  of  computing  the  two-dimensional  DFT  is  to  use  linear  congruential  scan¬ 
ning  of  the  data.  The  two-dimensional  discrete  Fourier  transform  system  will  use  an  input  scanning  device 
and  a  one-dimensional  discrete  Fourier  transform  device  as  shown  in  Figure  1 3.  The  two-dimensional  trans¬ 
form  block  size  N  j  by  N2  is  chosen  such  that  N  j  and  N  ->  are  relatively  prime  integers  (i.e.,  they  must  have 
no  common  divisor).  The  one-dimensional  Fourier  transform  device  has  a  block  length  of  N  =  N  |  N->.  The 
purpose  of  the  input  scanning  device  is  to  so  order  the  input  data  that  the  one-dimensional  Fourier  trans¬ 
form  of  the  length  N|Ns  serial  data  string  is  identical  to  an  N|  by  N-»  two-dimensional  Fourier  transform 
of  the  N  |  by  N2  input  data  samples.  If  desired,  an  output  scanning  device  may  also  be  used  to  provide  the 
transform  output  points  in  normal  order.  The  required  scan  may  be  derived  from  the  representation  115) 
of  a  one-dimensional  discrete  Fourier  transform  matrix  as  a  direct  product  matrix.  The  one-dimensional 
DFT  in  equation  (21)  is  equivalent  to  the  two-dimensional  DFT  in  equation  (22)  when  N  =  N|Ns  and  Nj. 
Ns  are  relatively  prime 
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for  n.  k  =  0.  I 


(21) 


Output 


Figure  13.  Two-dimensional  Discrete  Fourier  Transform  Device 
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forn^.kj  =0,  1.  ...  Nj  -  l  and  n-»k-,  =  0.  1.  .  .  .  Ns  -  1 

In  order  to  make  the  twolransforms  equivalent,  it  is  necessary  to  find  a  pair  of  one-to-one  functions 
n(n  | ,  n->)  and  k(k  ( ,  ki)  such  that 


n  j  k  |  n->k-> 

N,  +  n"T 


-i.  nik.  n-ik-, 

-TT-^-TT-1  +  -T 7^  (Modulo  1) 

N  N  |  Ns 

(231 

or 

nk  =  n|k|N->  +  n-,k->N|  (ModuloNjNj) 

(24) 

This  may  be  accomplished  by  letting 

n(n  j.  n->)  =  njN-i  +  n-»N  j  (Modulo  N) 

(25) 

k(k|.  k2)  =  k|U|N2  +  k^UiNj  (ModuIoN) 

(26) 

where  the  constants  Uj  and  U->  are  the  solutions  of 

N-)Uj  =  1  (Modulo  N] ) 

(27) 

N  |  Us  =  1  (Modulo  N-») 

(28) 

107 


I  quations  i  2  ^  >  anti  (  2h)  will  have  solutions  if  anti  only  it  N  j  anti  N  ^  arc  mutually  prime.  I  ndor  this  coiidc 
lion,  the  mappings  described  by  equations  (  25 )  and  (  2b)  will  sal ist >  the  requirement  o!  etiiuiion  i  25)  I  he 
linear  congruential  scan  prescribed  hy  I  qua  lion  (  25)  may  also  be  used  to  perform  two-dimensional  convo¬ 
lution  or  crosscorrelation  using  an  ordinary  one-dimensional  transversal  filter  or  crosscoi relator 


TWO-DIMENSIONAL  CZT 

Once  the  data  is  in  a  two-dimensional  format,  with  simultaneous  serial  access  to  all  the  rows,  it  may 
be  transformed  in  the  “horizontal”  direction  by  the  structure  shown  in  Figure  14  The  individual  one¬ 
dimensional  discrete  chirp  filters  and  discrete  chirp  generator  of  Figure  14  may  be  implemented  using  CCDs 
or  hybrid  digital  correlators.  An  acoustic  surface  wave  device  with  multiple  input  taps  may  he  used  to 
access  a  column  of  the  partially  transformed  output  in  a  single  shift  time  of  the  partial  transform  device. 
With  appropriate  coding  of  the  surface  wave  column  access  device,  it  may  also  perform  the  discrete  chirp 
premultiplication  and  the  discrete  chirp  convolution  of  a  DFT  in  the  “vertical"  direction.  A  complete  two- 
dimensional  CZT  architecture  is  shown  in  Figure  15.  and  the  required  coding  for  the  column  access  wave- 
device  is  shown  in  Figure  lb.  The  complex  arithmetic  may  be  implemented  as  described  previously  .  A 
balanced  mixer  may  be  used  for  the  fast  multiplier  required  for  the  vertical  transform.  Lower  speed  variable 
transconductance  multipliers  may  be  used  in  the  horizontal  partial  transform. 


Figurr  14.  Two-dimensional  Partial  Chirp-Z  Transform 
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Figure  15.  Hybrid  Implementation  of 'Two-dimensional  CZT. 
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Figure  16.  Tap  Weights  and  Structure  for  Acoustic  Surface  Wave  Combined 
Demultiplexer.  Chirp  Multiplier,  and  Discrete  Chirp  Filter 
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MODULAR  CZT 


Two  methods  for  combining  N  •»  chirp-7,  transform  (CZT)  modules  of  length  N  |  to  perform  a  dis¬ 
crete  Fourier  transform  (DFT)  of  length  NjNs  are  described.  The  first  method  uses  an  auxiliary  parallel- 
input.  parallel-output  DFT  device  of  si/.e  N  ->  and  allows  the  transform  of  size  N  | Ns  to  be  performed  in  the 
same  time  required  for  a  single  CZT  module  to  perform  a  si/e  N  |  transform.  The  second  method  uses  an 
auxiliary  parallel-input,  serial-output  DFT  device  of  si/.e  N  s.  If  the  second  method  is  implemented  entirely 
in  a  single  technology,  such  as  with  CCDs,  it  performs  the  si/e  N  |  N  s  transform  in  N  s  times  the  amount  of 
time  required  for  a  single  CZT  module  to  perform  a  size  N  |  transform:  if  N  s  is  a  composite  number,  say 
Ns  =  M  |  Ms.  the  second  method  also  permits  the  same  hardware  to  perform  M  |  simultaneous  transforms  of 
length  N | Ms. 

A  one-dimensional  discrete  Fourier  transform  may  be  written  as  a  partial  transform  of  a  doubly  sub¬ 
scripted  representation  of  the  data,  followed  by  a  pointwise  multiplication,  followed  by  a  second  partial 
transform  as  shown  in  equations  (29)  -  (34): 
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Figures  17  and  17a  show  modular  CZT  implementations  which  follow  from  equation  ( 34).  The  individual 
CZ7  subsystems  shown  in  Figures  I  7  and  1 7A  would  be  similar  to  tire  CZT  implementations  previously 
described.  The  parallel  DFT  required  for  the  second  partial  transform  in  Figure  17  may  be  implemented  as 
combination  of  summers  and  attenuators.  This  is  shown  in  Figure  IK  for  Ns  =  2.  In  general,  the  attenua¬ 
tion  factors  are  complex.  A  complete  double  length  CZT  is  shown  in  Figure  19.  Unfortunately,  a  parallel 
DF  f  implementation  of  this  type  becomes  unwieldly  if  the  dimension  Ns  is  very  large. 
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Figure  17.  Organization  ot  Modular  CZT 
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Figure  I  "A  Alternate  Organization  ot  Modular  CZT 
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I  It  cm:  equations  It  jvc  hccn  used  l<>  design  a  transform  ilex  tee  in  whiJi  signals  jic  'Ini  led  l  In  on  eh  t wo  <  <  I  > 
delav  line'  .it  d  if tci cn I  speeds  Alternatively  .  it  the  factors  in  equation  (  are  nitcipicted  j'  two  w.o.  - 
propagating  in  opposite  directions  relative  to  the  function  to  he  transformed  it  mav  he  seen  that  t Ik  o  i 
tare  of  I  igurc  2(i  also  perlornis  a  diserete  1  -miner  transform  with  speed  ..ompaiahlc  to  that  ot  a  (  /  I  \ 
surlaee  wave  device  module  which  implements  the  triple  product  convolutional  expiation  (d't  lias  heen 
built  hv  Reeder  |  l(>|  at  United  Aircraft.  A  schematic  is  shown  in  I  mure  21 


TWO-DIM  tN  SION  A  L  IK  T 

A  two-dimensional  IX  I  may  he  computed  as  a  two-dimensional  Fourier  transform  ot  a  symmetri¬ 
cally  extended  data  block  In  order  to  minimize  the  required  tiller  length,  the  scanning  ot  the  data  block 
which  has  been  symmetrized  for  the  two-dimensional  OIXT  will  now  be  considered  in  detail  I  he  data 
block  i .  extended  using  the  double  mirror  symmetry  defined  in  equation  (.US) 


g<  *  n  | .  t  n  s )  =  g(  n  | ,  ns) 


l  .Us  i 


The  two-dimensional  DPT  of  the  extended  data  block  defines  the  two-dimensional  DC  T  of  the  original  data 
block  given  in  equation  ( .M).  with  M  j  =  2N  j  -  I ,  M  •>  =  2N  ■>  -  1 
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figure  20  Parallel-Input,  Serial-Output  OZT  Using  Multi-Port  Transversal  Filter 
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Figure  21  Diode-Correlator  DFT  Modulo 

Since  the  indices  are  onl\  defined  Modulo  M  |  and  M  ->  respectively,  the  summation  limits  in  (  3‘)  i  are  re. ills 
no  different  front  those  in  (22)  If  the  scan  is  defined  hy  equations  (  25 )  and  (  2(i).  then  the  s\  niinetn  ol 
the  corresponding  one-dimensional  sequence  is  shown  in  equation  <  40). 

f_n  =  g(  -n  | .  -n  i )  =  g(  n  | .  n  i  =  t  n  1 40 1 

The  toroidal  scan  of  the  extended  data  block,  which  can  be  obtained  In  repeatedly  scanning  points  of  the 
original  data  block,  is  illustrated  in  Table  1  for  a  block  si/e  of  2  by  3.  The  numbers  in  the  table  indicate  the 
scan  order,  while  the  letters  indicate  the  data  values.  Similar  results  have  been  obtained  for  a  mixed  OIK" 1 
by  EDCT  which  permit  the  use  of  a  square  block  si/e.  |  I  7  | 


SIMULTANEOUS  COMPUTATION  OF  THE  DFT  AND  THE  DCT 

The  close  relationship  between  the  DFT  and  the  DCT  permits  the  use  of  common  modules  to 
simultaneously  compute  both  transforms.  This  may  be  accomplished  most  simply  when  an  EDCT  is  com¬ 
puted  using  DFT  modules.  The  sum  in  the  EDCT  defining  equation  (  lHa)  may  be  interpreted  as  a  length 
2N  DFT  of  the  extension  of  the  function  g  by  N  zeros.  This  leads  to  the  configuration  shown  in  Figure 
22.  Alternatively,  if  the  odd  and  even  frequencies  in  the  /.cro-filled  DFT  are  considered  separately,  they 
may  be  computed  using  length  N  DFT  modules  as  shown  in  Figure  23. 
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Figure  22.  Computation  of  the  DFT  and  EDCT  Using  a  Single  Length  2N  DF1  Module 
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Figure  23.  Computation  of  the  DFT  and  EDCT  Using  Two  Length  N  DFT  Modules 
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CONCLUSION 


It  lias  been  shown  tli.it  .1  small  inmibei  of  multipurpose modules  can  pioudc  .1  highly  patullcl 
structure  with  nimimal  control  overhead  tor  the  coniput.ition  o|  many  linear  Hailstorms  I  lie  tr.msloi ms 
include  one-  and  two-dimensional  discrete  Fourier  translornis  and  discrete  cosine  trails!  orms  I  he  b.isu 
modules  may  lie  chosen  to  he  multipliers,  serial  access  memories,  and  discrete  chirp  tillers  01  ciosscoirelaloi  s 
Current  acoustic  surface  wa\e  and  (CL)  technologies  permit  small,  low  power,  lightweight,  high  speed 
implementations  of  the  required  modules  and  permit  real-time  solutions  to  highly  demandin'.’  signal  pov 
essmg  tasks  including  real-time  video  data  compression.  |  I  <S  |  spread  spectrum  conununii  at  1011  and  high 
resolution  radar  signal  processing. 
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APPENDIX  H 


INVESTIGATING  ADAPTIVE  TRANSFORMS  IN 
INFORMATION  PROCESSING  APPLICATIONS 


PREFACE 


This  report  describes  work  performed  during  the  five  month 
period  from  June  1,  1975  through  November  27,  1975.  The  investi¬ 
gation  was  funded  by  the  Naval  Undersea  Center,  San  Diego,  Calif¬ 
ornia  under  contract  no.  N  66001-75-0-226  MJE,  and  was  conducted 
at  UCLA  under  the  direction  of  Judea  Pearl  as  Principal  Investi¬ 
gator.  The  team  engaged  in  this  study  consisted  of:  Judea  Pearl, 
Massih  Hamidi  and  Yechiam  Yemini. 

The  continuous  guidance  and  encouragement  of  Harper  Whitehouse, 
Jeff  Speiser  and  Robert  Means  of  the  Naval  Undersea  Center  deserve 
the  major  credit  for  the  accoinpl  i  shments  reported  here. 


ABSTRACT 


The  performances  of  the  Discrete  Cosine  Transform  (CCT)  and 
the  Discrete  Fourier  Transform  ( DFT )  were  analyzed  and  compared 
in  signal  processing  applications.  Conditions  for  the  asymptotic 
optimality  of  the  DFT  and  DCT  were  established.  The  superiority 
of  the  DCT  over  the  DFT  was  established  for  Markov-1  signals.  A 
tool  for  studying  asymptotic  behavior  of  transform  was  developed 
using  numerical  quadrature  analysis. 
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SUMMARY  OF  WORK  PERFORMED 


Guided  by  Naval  Undersea  Center  directives  the  emphasis  of 
our  investigation  has  been  to  analyze  and  compare  the  performance 
of  two  practical  transform  techniques,  the  Fourier  (DFT)  and  the 
Cosine  (OCT)  transformation,  in  image  processing  applications. 

Our  study  has  resulted  in  several  new  results  with  both  theoreti¬ 
cal  and  practical  significance. 

1.  While  it  is  well  known  that  the  Fourier  coefficient  of  sta¬ 
tionary  finite-duration  continuous  signals  are  asymptotically 
uncorrelated,  correlation  properties  of  finite-dimensional, 
discrete  Fourier  Transform  remained  an  open  question.  We  have 
established  the  following  facts: 

a.  If  the  covariance  sequence  is  summable  the  magnitude  of 
every  off-diagonal  covariance  element  converges  to  zero 
as  N  +  ». 

b.  If  the  covariance  sequence  is  only  square-summabl o  the 
magnitude  of  the  covariance  elements  sufficiently  far 
from  the  diagonal  converges  to  zero  as  N  -*  °°. 

c.  If  the  covariance  sequence  is  square-summabl e  the  weak 
norm  of  the  matrix  containing  only  the  off-diagonal  ele¬ 
ments  converges  to  zero  as  N  -*■  ®. 

d.  If  the  covariance  sequence  is  summable  the  weak  norm  of 
the  matrix  containing  only  the  off-diagonal  elements  con¬ 
verges  to  zero  at  least  as  fast  as  -Jp  . 


2.  It  was  conjectured  that  the  performance  of  the  Cosine  trans¬ 
form  (DCT)  is  superior  to  that  of  the  Fourier  transform  (DFT). 
The  fact  that  DCT  was  found  more  compatible  with  the  hard w are 
configuration  of  the  Image  Processing  group  at  NUC  called  for 
an  analytical  examination  of  this  conjecture.  He  were  success¬ 
ful  in  establishing  the  following  results: 

a.  The  DCT  is  asymptotically  equivalent  to  the  Karhunen-Lo£ve 
transform  (KLT)  of  Markov-1  signals  and  the  rate  of  conver¬ 
gence  is  (similar  to  the  DFT)  on  the  order  of  N"1/2. 

b.  The  DCT  offers  a  better  approximation  to  the  KLT  of  Markov- 
1  signals  than  the  DFT  for  all  values  of  N  and  p. 

c.  The  DCT  is  asymptotically  equivalent  to  the  KLT  of  all 
finite-order  Markov  signals. 

d.  The  analysis  of  asymptotic  properties  of  discrete  trans¬ 
form  can  be  simplified  substantially  using  numerical 
quadrature  techniques. 

A  detailed  description  of  these  results  are  containea  in  the 
following  three  appendices. 
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ON  THE  RESIDUAL  CORRELATION  OF  FINITE¬ 
DIMENSIONAL  DISCRETE-FOURIER- 
TRANSFORMS  OF  STATIONARY  SIGNALS 

Massih  Hamid i 
and 

Judea  Pearl 


Published  in  IEEE  Transactions  on 
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ABSTRACT 


The  covariance  matrix  of  the  Fourier  coefficients  of  N-sampled 
stationary  random  signals  is  studied.  Three  theorems  are  established. 

1.  If  the  covariance  sequence  is  sunmable  the  magnitude  of  every 
off-diagonal  covariance  element  converges  to  zero  as  N  +  ». 

2.  If  the  covariance  sequence  is  only  square-sutrmable  the  magnitude 
of  the  covariance  elements  sufficiently  far  from  the  diagonal 
converges  to  zero  as  N  -*■  «*. 

3.  If  the  covariance  sequence  is  square-summable  the  weak  norm  of 
the  matrix  containing  only  the  off-diagonal  elements  converges 
to  zero  as  N  -*  <*>. 

The  rates  of  convergence  are  also  determined  when  the  covariance 
sequence  satisfies  additional  conditions. 
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1.  INTRODUCTION 

It  is  well  known  that  the*  Fourier  coefficients  of  stationary  finite-dura¬ 
tion  continuous  signals  are  asymptotical  ly  uncorrolated.  Root  and  Pitcher^  ^ 
have  shown  that  the  cross-correlation  between  any  two  Fourier  coefficients  ap¬ 
proaches  zero  as  the  signal  duration  grows  infinite.  Consequently,  problems 
involving  stationary  stochastic  processes  are  often  treated  by  approximating 
the  original  processes  by  Fourier  series  with  uncorrelated  random  coefficients. 

In  applications  involving  digital  signal  processing  the  continuous  signals 
are  sampled  at  a  finite  number  N  of  equally  spaced  points  in  time  and  are 

treated  as  N-dimensional  vectors.  Root  and  Pitcher's  results  coupled  with  the 

C21 

computational  convenience  of  the  Fast-Fourier-Transfornr  J  gave  rise  to  a  num¬ 
ber  of  applications  whereby  the  finite  Fourier  transform  of  the  signal  vector 
is  taken  and  its  components  are  treated  as  uncorrelated  random  variables.  For 
example,  in  digital  transmission  of  pictures  and  voice  it  is  a  common  practice 

to  assign  each  Fourier  component  a  digital  code  which  is  independent  on  the 

T31 

magnitude  of  the  other  Fourier  coefficientsL  . 

The  Fourier  coefficients  of  a  signal  N-vector  x^  =  (xQ,  x^,  ...  xN1)T 
are  uncorrelated  if  and  only  if  the  covariance  matrix  of  x^  is  a  circulant 
matrix^;  i.e.,  if  £(x-x.)  is  a  function  only  of  (i-j)  mod  N  .  Since  such 

'  J 

circular  symmetry  is  very  rare  in  actual  processes  the  covariance  matrix  of 
the  Fourier  coefficients  will  contain  off-diagonal  elements  whose  magnitudes 
affect  system  performance. 

It  is  generally  believed  that  the  magnitude  of  each  off-diagonal  element 
and  their  cumulative  effect  both  converge  to  zero  as  N  -+  »  .  In  this  paper  we 
establish  conditions  under  which  these  convergences  take  plac’  and  derive  ex¬ 
pressions  for  the  rate  at  which  the  residual  cross-correlation  decays  with  N  . 


2.  THEOREMS  ANO  PROOFS 

T 

Consider  a  sampled  sequence  x^  =  [x^,  x-j,  -  xfJ  j  ]  of  a  stationary 

stochastic  process  with  a  Tosplitz  correlation  matrix  Ejx^,  xl}  satisfy¬ 

ing 

(V..  =  TN(  I  >  -  J  I  )  =  t(  |  i  -j  | )  i  >  j  -  0  ,  1  ,  - N-l 

^  J 

where  x  designates  the  complex  conjugate  of  x  . 

The  discrete-Fourier-transform  (0  FT)  of  is  a  vector  -  [y^ ,  ••••  -j  UT 
defined  by  y^  =  FN  xN 

where 

[Fn]  =  N'V2  W™  m,n  =0,  1,  *•••  N-l 

m,n 

and 

«N  -  e,2’/N 

We  study  the  behavior  of  the  off-diagonal  elements  of  the  matrix 

^  =  FN  tn  fn 

for  large  N  ,  as  well  as  the  norm  of  where 

CN  a  diag.j  E(|yQ|2)  ,  E(ly-,!2)  ■•••  E(|yrM!2) 

The  motivation  for  studying  this  norm  lies  in  the  fact  that  in  many  signal  pro¬ 
cessing  applications  the  performance  degradation  caused  by  the  residual  corre¬ 
lation  can  be  upper  bounded  by  monotonic  decreasing  functions  of  | T^  -  C^j^. 

Theorem  1 

If  t(*)  is  summable  then 

E{ym  yn)  =  0  m  t  n 

N  -»•  ® 

Moreover,  lim  E{|y  |2}  =  2  \  t{u)  -  t(o)  =  \  t(u) 

N  -  -  m  u,0  u*— 
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Proof:  By  definition 


.  N-l  N-l  r  ,  N-l  ,  v  t 

MVnl-  »■'  i  r  T„  «inr  ■  f’  i  I  < 

*,=o  r=o  1  nr  i.-o  u-t-(N-l ) 


"Ut((u| 


In  order  to  obtain  an  expression  with  only  single  summation  we  use  Abel's  Trans¬ 
formation  Formula  for  Partial  Summation^: 

l  ^  ’  jmAk(bk  -  *W  +  flnb„  '  \.-lbm 

where 


A 


k  "  E  aj 

k  j-o  J 


for  k  ^  0  and  A  ^  =  0 


Letting  a£  =  W^m"r'^ 


and  b  =  l  wlju  t(u)  we  obtain 

U=£-(N-1 )  N 


NE(v„t  ■  l 

1=0 


N-l  M 


i  wi;'™-1'1 


k=o 


Wj*  [t(N-t)  -  t(t)] 


for  m  f  n  ,  and: 


N-l 


N-l 


,rrur 

’m1  J  ^  'N 

£  =  0  U=0 


N  E{|ymr}  =  l  CWN-O-tWl  +  N  l  w™  t(u) 

From  equation  ( 1 ) : 


:!v„> 


<  N 


-1 


N-l  ( £.-1 


l 

1=0 


l  H 


k(m-n) 


k=o 


The  first  term  can  be  bounded  by  N"^  l  fc|tU)| 

£=o 

For  the  second  term  we  use  (for  m  f  n) 

V  =  o  =  N  I  1  +  V 

k=o  N  k=o  N  k=N-u  N 


N 

N-l 


Wjj*  t(z) 


+  N 


-1 


N-l ( 2.-1 

I 


0) 


(2) 


2=0 


k=o 


and  obtain 


,-l 


N-l  f 2.-1  , 

l  l  «n 

£=olk=0 


M[Jlt(N-£) 


.-1 


N 


<  N"  l  u  1 1 ( u )  | 
u=l 
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Hence: 


E{yn,yn}|  i  2N']  l  *|tU)l 

^  0 


Since  t(t)  is  summable;  given  e  >  0  ,  there  exists  an  integer  M  such  that 


l  |t(Of  <  | 

*=M  ^ 


Thus: 


|E<Wn>|  i 


,  M-l  N 

N"1  l  *|t(*)l  +  I  | t(£) | 

£=o  £=M 

-1  M'1  e 

Now  choosing  N  large  enough  to  have  N  [  t|tU)l  j  establishes  the 

1=0 

first  part  of  our  proof. 


From  equation  (2),  we  obtain,  after  some  computation 


N-l 


N 


E{|yJ^}  -  l  W™  t(u)  +  l  W-mut(u)  -  N"1  l  u 


N-l  r 


'm1 


N 


u=o  "  u=l 

taking  the  limit  of  each  term,  we  have: 

I  N-l 


u=l 


W™  +  W-mu|  t(u) 


1  im  N 


-1 


N 


l  uWjMUt(u) 
u=o 


N-l 


<  lim  N_1  l  u|t(u)|  =  0 

N  -*•  »  u=o 


also,  it  can  be  shown  that 
N  jtmu 


lim  l  Vfcj  u  t(u)  =  l  t(u) 
N  ->  co  u=o  u=o 


(3) 


Hence: 


lim  E{ |ym| =  2  l  t(u)  -  t(0)  =  l  t(u) 


N  ® 


u=o 


u=-° 


Q.E.D. 
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To  show  that  (3)  is  valid  notice  that  t(u)  being  suinnable.  given  any  r  >  0 
there  exists  an  integer  K  such  that 


V  t ( u)  <  |  and  l  j  t ( u)  |  W 
u=K  c  u-K 


l  t(u)  -  lim  l  W±mu  t(u) 
u=o  N  ■>  »  u=o 


oo  K-l 

N  +m„ 

< 

I  t(u)  -  l  t(u) 

+ 

1im  I  WN  t(u) 

u=o  u=o 

N  ->■  «  u=K 

CO  N 

l  t(u)  +  lim  l  1 1 ( u) |  <  e 

u=K  N  -*■  ®  u=K 


Theorem  2 

If  t(t)  is  square  summable,  lim  E{y_yn}  =  0  for  all  elements  such 

N  -v  o°  m  n 

that  N  I m-n |  >  e  >  0  . 


Proof: 


From  equation  (1 ) : 


EfV„! 


1  -  W 


(m-n) 


-1 


N-l  , 

y  wni  - 
An  un 
2=0 


[t(N-2)  -  t(0] 


Applying  the  Cauchy-Schwarz  inequality: 


Etv„( 


<  N 


-1 


1  -  W 


(m-n) 


<  2/2  W 


N 

1/2 


-IfN-l 

I 

k2=0 


un£  ,.nu 
WN  "  WN 


2'1/2fN-l 


i  it(N-£)  -  t(or 
£=0 


1/2 


1  -  e 


i  (m-n)2Tr/N 


-1 


'  N 

l  |tU)l‘ 
1=0 


1/2 
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1 


Hence,  whenever 


N"1 | m-n | 


>  e  >  0 


Theorem  3 


lim  I 

N  ->  cJ 


summable. 


t  i  N-l  _  ? 

"  mi J£{vJ 

m,n-o 
m  1  n 


1/2 


0 

Q.E.D. 


0  ,  if  t ( u)  is  square 


Proof: 

r  5 1 

From  equation  (1)  we  obtain  with  some  computations1  J 

,  N-l  o  N-l  .  P 

N2  l  IMy^n!  =  N  l  (Hz  -  £2)|t(N-£)  -  t(£)|2 
m,n=o  ?=o 

m  f  n 

Hence 

lTN  '  SI2  =  N_1  V  !E^Vn^  =  N"2  V(N^  -  £2)  I t(N-pJ  -  tU)!2 

‘  '  flTjft-O  £=0 

m  f  n 

<  4N"2  V  (Mu.  -  £2)|t(£)|2 
£=0 


Let  e  >  0  ,  arbitrary,  be  given.  Since  t(£)  is  square  summable  there  exists 

00  p 

an  integer  P  such  that  7  |t(si)|  <  §•  .  Thus 

a-p  0 


TN  -  CK  i  4 


N'2  Y  (Nt  -  J2)|t(jt)|2  +  l  |tU)|2 
£=0  £=P 


p  P-1 

Choosing  N  large  enough  to  have  N  £  (Nt 

£=0 


i2)|tU)|2 


<  |  ,  yields 


TN  '  CN 


<  E 


Q.E.D. 


The  convergence  rates  of  both  the  off-diagonal  elements  and  the  norm  of 
Tn  -  can  be  obtained  under  certain  conditions  from  the  proofs  of  theorems 
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1-3.  These  will  he  listed  in  the  following  three  corollaries: 


Corollary  1  -  If  ct(».)  is  summable  then  Efy^^y  }  =  o(N)  for  m  1  n  . 

Corollary  2  -  If  t.(e)  is  square  summable  then  =  o(N^)  for 

all  | m-n j  =  0(H) 

Corollary  3  -  If  i}^  t(n)  is  square  summable  then  |T^  -  C^j  =  o(N^2)  . 

Note  that  these  conditions  on  t[&)  are  satisfied  for  most  processes  en¬ 
countered  in  practice,  e.g.,  finite  order  Markoff  or  moving-average  processes. 

3.  CONCLUSIONS 

The  theorems  and  corollaries  established  in  this  paper  constitute  a  gen¬ 
eralization  of  Root  and  Pitcher's  results  to  the  case  of  discrete-time  signals 
and  provide  guidelines  for  selecting  the  proper  block  length  N  in  Fourier- 
signal  -processing  applications. 
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APPENDIX  II 


COMPARISON  OF  THE  COSINE  AND  FOURIER 
TRANSFORMS  OF  MARKOV-1  SIGNALS 

Massih  Hamidi 
and 

Judea  Pearl 


Submitted  for  publication  in  IEEE  Transactions  on 
Acoustics,  Speech,  and  Signal  Processing 


I.  INTRODUCTION 


In  a  recent  paper  Ahmed,  et  al.^  propose  a  new  transform  called  Discrete 
Cosine  Transform  (DCT)  and,  based  on  empirical  evidence,  conjecture  that  its  per¬ 
formance  is  closer  to  the  optimal  Karhunen-Loeve  transform  (KLT)  than  the  other 
commonly  used  transforms  (i.e..  Discrete  Fourier,  Walsh-Hadamard,  Haar).  Means, 
et  al.^  actually  used  the  DCT  for  encoding  TV  pictures  in  real  time. 

Pearl  showed^  that  for  a  signal  statistic  characterized  by  a  covariance 
o 

matrix  T  ,  | T  -  Ty |  (to  be  defined  later)  constitutes  a  measure  of  performance 

for  a  transform  U  ,  in  the  sense  that  the  error  bounds  (in  coding  and  filtering) 

2 

are  increasing  functions  of  |T  -  Ty [ 

The  purpose  cf  this  investigation  is  to  determine  the  relation  between 
|T  _  t  | c  (the  norm  obtained  using  the  discrete  cosine  transform)  and  |T  -  T p ( 
(the  norm  obtained  using  the  discrete  Fourier  transform  (CFT)),  thus  testing  the 
conjecture  of  Ahmed  and  his  collaborators. 


II.  DEFINITIONS  AND  NOMENCLATURE 

I 

Let  T  be  a  Toeplitz  matrix  and  U  an  orthogonal  transform.  Let  T 

UTtT^  be  the  representation  of  T  in  the  new  basis,  and  T^  =  diag 

(T11’  T22’  "*»  Tii’  •**»  V  '  We  define  Tu  t0  be  the  ^presentation  of 

T.  in  the  first  basis,  i.e., 
li 

T(J  -  Ty  u 


and  |T  -  Ty |  the  Hilbert-Schmidt  norm  of  T  -  Ty  ,  i.e., 

M-l 


lT  "  V 


1_ 

M 


KT-  V 


(m,n=o 


mn 


The  cosine  transform  representation  of  a  Toeplitz  matrix  T  is  given  by 
CTC"^  where  C  is  a  MxM  matrix  defined  by 
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"OJ 


2  1 

H/7 


ckj  =  ft  cos  M  (2J+]) 


A  simple  algebraic  manipulation  shows: 


*  *  o 

CC  =  C  C  =  |  I 


j  =  0,  1,  •••,  M-l 

k  =  1,  2,  ,  M-l 

J  =  0,  1,  M-l 


(where  X  indicates  the  complex  conjugate  transpose  of  X  ).  Hen  ■»: 


c-'  -  fc* 


In  contrast,  the  DFT  is  defined  by  a  unitary  matrix  F  where 

Fkj  =  M"1/2  exp|i  kjj  k,j  =  0,  1 ,  •  •  • ,  M-l 


III.  COMPARISON  OF  |T  -  Tj2  AND  |T  -  Tp | 2 

For  any  orthogonal  matrix  U  (e.g.  U  =  C  or  U  =  F  )  we  have 

M-l 

l 

m=o 

2 


IT  -  T.  .1  c  =  |  T*  -  Ty  |  2  =  |t'|2  -  1  "j'  KUTU-1) 


'U' 


mm1 


M-l 


■  Pl2-J  I 


m=o 


i.e.:  The  higher  the  norm  of  the  diagonal  vector  of  the  transformed  matrix, 

the  lower  |T  -  Tyj  and  the  better  the  transform.  Hence,  to  compare  |T  -  Ty|‘ 
and  |T  -  Tp|  it  suffices  to  compare 


M-1  .  -  M-l  ,  ? 

l  l(CTC  >J  and  l  KfTF  )J 

m=o  m=o  1 


We  consider  matrices  of  the  form 
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T  = 


1  P  p‘ 

P  1 


M-1 


M-2 


M-1 


which  represent  covariance  matrices  of  Markov-1  signals, with  0  _<  p  _<  1 
being  the  covariance  coefficient  between  adjacent  samples. 

Clearly,  for  p  =  0  and  p  =  1  the  cosine  and  Fourier  transforms  are 
equivalent,  since  T  is  diagonal  in  both  representations.  For  an  intermediate 
value  of  p  we  obtained: 


and 


( CTC  1 ) 

'oo 


1  +  P  2  p(1  -  P^) 
1  -p~M  (W? 


(CTC'1  ) 


mm 


1 

2 


pp  -  Htf). 

M 


for  m  f  0,  where  a  =  hitt/viM. 

An  elementary  (but  tedious)  computation,  leads  to: 


M-1  , 

l  | (CTC-1) 
m=o 


M(l+p2)  4p2 

,  2  "  77 ~T? 

1  -  p  (1-P  ) 


ry  f  *1  2M\  2  n 

•  [3(1+P2)  +  4P] 

Md-p2r 


40;n-0")2 

M20-p)4 


for  M  =  2k  ,  k  >  1  (i.e.  M  even  ^4)  . 
Combined  with 

9  i.2  A  2#i  2M\ 

Mill2  =  m  1  P  -  .?P  il  -V-  1 
nm  m  2 

1  -  p  (1  -  p  ) 


we  finally  obtain  the  desired  norm: 
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«|T-TC| 


2 


2..2n  *  p2M1  .  2p2(1  -  p2") 

(I  -  p2)2  M(1  -  P2)3 


O  a  2 f 7  M \ 2 

(3  +  4p  +  3p2)  +  ^yl-~  -P-l 
M2(l  -p)4 


Note  that 

lim  | T-T  |  =  0 

M  •  >  oo 

implying  that  the  DCT  is  asymptotically  equivalent^  to  the  KIT  of  Markov-1 
processes.  Moreover,  since  for  large  M  and  p  j*  1  we  have 


| T-Tc |  =  /?  -^7  0(M"1/2) 

1-p 

we  conclude  that  the  degradation  in  performance  in  filtering  and  coding^  van¬ 
ishes  like  M ~^2  . 

In  order  to  calculate  |T  -  Tr|^  recall^  that  for  T..  =  t(|i-j|)  we 

*  I  J 

have 


(T  -  T  )  =  Jiii  [t(|H|>  -  t(M  -  I  HI)] 

ij 

and  substituting  t(|i-j|)  =  pi1'-'*!  we  obtain: 


M|T  -  Tp|2 


2p2(1  +  p2W)  2(1+p2)p2(1-p2M)  _  pM(M2-1) 
(1  -  p2)2  M(1  -  P2)3  3 


It  shows  that  the  asymptotic  behavior  of  |T  -  Tp |  for  large  M  is  identical 
to  that  of  | T  -  Tc |  .  Thus,  the  performance  difference  between  the  DCT  and 
the  DFT  must  vanish  like  .  Indeed,  for  large  M  one  obtains  the  posi¬ 
tive  difference 


IT  - 


IT 


ie? _ _ 

m2(i  -  P2)(i  +  p)2 


p  <  i 


indicating  that  the  cosine  transform  is  closer  to  optimal  than  the  Fourier 
transform  over  the  entire  range  of  0  <  p  <  1  . 
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For  moderate  values  of  M  we  should  examine  the  expressions  for  |T  -  Tp  | 
2 

and  |T  -  T^l  over  the  range  0  _<  p  _<  1  .  The  two  are  plotted,  in  a  normal- 
*  •  2 

ized  form,  in  Figure  1.  We  chose  J T  -  I [  as  a  common  normalizing  factor, 
where  I  is  the  identity  matrix,  and  so 

I T  -  I|2  =  — EM-1  “  MP2  +  P2M]  . 

M(l-p2)2 

It  measures  the  degree  of  cross-correlation  contained  in  the  unprocessed  sig¬ 
nal,  and  therefore  the  maximum  amount  of  decorrelation  that  can  be  accomplished 
by  any  transform  (i.e.  the  KLT).  The  ratio 


represents  the  fractional  correlation  left  'undone1  by  a  transformation  U  . 

Figure  1  shows  that  for  M  =  8,  16,  64,  and  for  the  entire  range  of 
2  2 

0  <  p  <  1  ,  [T  -  Tp|  is  higher  than  f  T  -  Tcf  .  The  difference  between 
the  two  are  quite  noticeable,  occasionally  reading  a  ratio  of  2  :  1. 

CONCLUSIONS 

We  established  that  the  DCT  is  asymptotically  equivalent  to  the  KLT  of 
Markov-1  signals  and  demonstrated  that  the  rate  of  convergence  is  on  the  order 
of  M~^2  .  |  T  -  Tj-|  2  is  shown  to  be  smaller  than  |T  -  Tp|  2  for  all  values 
of  M  and  ,  i.e.,  the  Discrete  Cosine  Transform  offers  a  better  approxima¬ 

tion  to  the  Karhunen-Lbeve  transform  of  Markov-1  signals  than  the  Discrete 
Fourier  Transform. 
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I.  INTRODUCTION 


Discrete  unitary  transforms  are  used  extensively  in  digital 
signal  processing.  Both  voice  and  images  were  successfully  en¬ 
coded  and  filtered  using  Fourier,  Cosine,  Walsh,  Haar,  and  Kar- 
hunen-LoSve  expansions^  All  unitary  transformations  are 

information  preserving  and  no  bandwidth  reduction  results  from 
the  application  of  the  transform  to  the  signal.  Instead,  its 
beneficial  effect  lies  in  redistributing  the  variance  associated 
with  each  signal  sample  into  almost  uncorrelated  variables,  thus 
permitting  the  transform  coefficients  to  be  processed  (e.g.,  co¬ 
ded  or  filtered)  individually^^.  When  the  signal  statistic  is 
known  one  can  find  an  optimum  transformation ,  the  Karhunen-Logve 
transform  (KLT),  which  totally  decorrelates  the  transform  coeffi¬ 
cients  and  results  in  an  optimal  performance. 

Real-time  implementation  of  the  Karhunen-Lo&ve  expansion  suf¬ 
fers  from  three  major  drawbacks:  1)  A  large  amount  of  sampled- 
data  and  computational  effort  is  required  for  estimating  the  di¬ 
rection  cosines  of  the  Karhunen-LoSve  basis  vectors.  2)  The  task 
of  projecting  an  N-component  incoming  signal  vector  into  its  prin- 
cipal  coordinates  requires  N  computer  multiplications.  3)  A  new 
NxN  matrix  must  be  calculated  for  each  change  in  the  statistical 
properties  of  the  environment. 

The  advantages  of  other  transform  techniques  (e.g.,  Fourier, 
Cosine,  Walsh)  lies  in  the  fact  that  they  possess  fast  computation 
algorithms  and  that  they  employ  a  fixed  set  of  unitary  matrices 
independent  of  the  signal  statistics.  Their  performances,  however 


should  be  judged  by  the  degree  to  which  each  transform  approxi¬ 
mates  the  Knrhunen-Lodve  expansion  in  decorrelating  the  signals. 

A  good  measure  of  the  degree  of  correlation  still  remaining 
after  the  application  of  a  specific  transform  is  given  by  the 
norm  of  the  matrix  containing  the  off-diagonal  covariance  ele¬ 
ments  of  the  transformed  coefficients.  This  norm  was  shown^-^ 
to  control  the  performance  degradation  resulting  from  residual 
correlation  in  both  coding  and  filtering.  When  this  norm  is 
small  one  is  assured  that  the  performance  degradation  will  re¬ 
main  below  a  tolerated  level  and  therefore,  the  behavior  of  this 
norm  for  large  N  governs  the  proper  selection  of  the  signal  di¬ 
mension  N .  When  the  norm  generated  by  a  specific  transform  van¬ 
ishes  at  large  N  we  say  that  that  transform  is  asymptotically 
equivalent  to  the  KLT .  The  rate  at  which  the  norm  vanishes, 
however,  determines  quality  ranking  among  several  contending 
transforms . 


Consider  a  continuous  signal  which  is  sampled  at  N  time 
points  to  give  an  N-dimensional  vector  ,x.  Let  be  an  N-dimen- 
sional  unitary  transform  matrix  and  y_  A  Cf)  x  the  transformed 
signal.  If  x^  is  considered  a  random  vector  with  an  autocovari¬ 
ance  matrix  T  then  the  transformed  autocovariance  matrix  is 
*  * 

CN  T  CN  (Cfj  is  the  adjoint  matrix  of  C^).  The  norm  determining 
the  degree  of  decorrelation  achieved  by  is  given  by 


2 

i  j 


Now  consider  a  sequence  of  such  signal  vectors  wit.i  increasing 
dimension  and  the  corresponding  sequence  of  autocovariance  matrices 


CO 

{Tm}  [N  will  be  called  the  block  size].  We  may  transform 
"  N=1 

each  T ^  using  a  transform  of  corresponding  dimension  taken 

CO 

from  a  transform  sequence  { C . , }  .  We  want  to  examine  the  de- 

N=  1 

gree  of  diagonal izat ion  of  the  transformed  sequence 
as  the  block  size  grows  to  infinity. 

If  the  signal  covariance  matrices  were  known,  it  would  be 
possible  to  compute  (numerically)  the  value  of  i- 
for  each  N,  observe  its  behavior  or  verify  asymptotic  equivalence. 
In  most  cases,  however,  a  system  of  transforms  must  be  decided 
upon  with  only  partial  knowledge  of  the  input  statistics.  Tech¬ 
niques  using  the  di screte-Fouri er  transform  (DFT),  or  discrete 
cosine  transform  (DCT)^,  for  instance,  are  often  expected  to 
process  stationary  signals  of  arbitrary  statistics.  Analytical 
techniques  must  be  employed  for  examining  asymptotic  behavior  of 
CN  Tn  CN  over  an  entire  class  of  covariance  matrices  .  A  basic 

* 

difficulty  impeding  such  analysis  is  that  each  element  of 
consists  of  sums  of  components,  where  both  the  component  values 
and  the  size  of  the  sum  are  changing  with  N.  The  expressions  re¬ 
ceived  for  those  sums  are  formidable,  which  make  it  almost  impos¬ 
sible  to  derive  analytic  results  regarding  the  asymptotic  behavior 
of  a  given  transform. 

This  paper  attempts  to  reduce  this  difficulty  be  deriving 
certain  limit  behavior  of  the  elements  of  the  transformed  matri¬ 
ces.  An  analytical  framework  is  developed  which  relates  the  prob¬ 
lem  at  hand  to  known  theories  of  analysis,  thus  making  it  possible 
to  utilize  known  limit  theorems  of  classical  analysis  to  derive 
limit  behavior  of  transforms. 


(CN  TN  Cn1n=1 
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2.  BASIC  DEFINITIONS 

Let  15^  be  the  space  of  symmetric  and  real  N-dimensional 

matrices.  On  we  use  two  norms,  the  weak  norm  and  the  strong 
H  N  - 1 

norm.  Let  AfS^,  be  a  matrix  with  eigenvalues  { >,  . }  £  .  The  weak 


norm  of  A  is  defined  by: 


9  •,  N-1  N-l  9 

\<  a  J  .z  i  I a i j i 

1=0  J=0  J 
norm  of  A  is: 


i  N-l 

N  ,  £JXi 


A  sup{  <Au,  u>  :  ||u||  =  1  u.CTRn} 


[where  <  ,  >  denotes  an  inner  product  on 


0  <  i  <  N  - 1 


Both  norms  are  invariant  under  unitary  transforms.  That  is,  if  C 


is  unitary  |A|  =  |C  A  C*|  and  ||A||  =  ||  C  A  C*j|.  Moreover,  ||A||N  >  |  A  |  ^ 
which  is  where  the  norms  acquire  their  names.  The  weak  norm  of  the 
off-diagonal  portion  of  the  transformed  covariance  matrix  (namely, 

| C  A  C*  -  DIAG(C  A  C*)| )  measures  the  degree  of  residual  correla¬ 
tion  between  the  transformed  signal  components  and  bounds  practi¬ 
cal  performance  degradation  resulting  from  residual  correlations^. 

To  be  able  to  consider  sequences  of  matrices,  such  as  the  au¬ 
tocovariance  matrices  of  signals  with  increasing  block  size,  we 
define  nets.  A  net  is  a  strongly  bounded  sequence  of  matrices 
{A,}”  »  such  that  for  every  N,  A^  is  a  NxN  matrix.  Two  nets. 


a  =  {An}  and  8  =  {Bn}  are  called 
lAN  '  bnL  T” °* 


equivalent  if 


A  matrix  class  A »  is  a  collection  of  nets.  The  N-section 
of  -A- is  the  collection  of  N-dimensional  matrices  wh;ch  belong  to 
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nets  in  A.  We  denote  the  N-section  of  A  by  >. 

Let  b  e  the  matrix  class  of  all  nets  of  positive  diagonal  ma¬ 

trices.  Ax  i  s  called  the  diagonal  class. 

Let  y  =  {C^}  be  a  net  of  unitary  transform  matrices  and  3  = 

{A^}  an  arbitrary  net.  Then  by  y  a  v  we  mean  the  transformed  net 

{CN  AN  ^  ^  1S  a  matr‘lx  class  we  denote  by  Y  A  y  the  trans- 

* 

formed  class.  y  jy y  we  call  the  y-spectral  representation  of  the 

★ 

class  A m .  y  >&  y  is  called  the  class  diagonal  in  y .  As  an  example, 
the  class  diagonal  in  the  net  of  Fourier  transforms  is  the  class  of 
circulant  matrices. 

Given  a  transform  net  y  we  consider  the  performance  of  y  on  a 
class  of  covariance  matrices  !f,  called  the  signal  class.  y  performs 
well  on  T,  in  the  sense  that  it  approximately  diagonalizes  the 
class  J,  if  every  net  y  x  y  in  y  Ty  is  asymptotically  equivalent 
to  a  diagonal  net.  This  motivates  the  following  definition;  if  A 
and  $  are  matrix  classes  A  is  said  to  be  an  asymptoti c  cover  of  ; 
if  for  any  net  8  €  (3  ,  there  is  a  net  a(A,  such  that  a  and  B  are 
asymptotically  equivalent.  We  use  the  notation  to  denote 

this  fact.  For  example,  the  class  of  circulant  matrices  is  known 
to  cover  the  class  of  Toeplitz  matrices^^.  Similarly,  the  class 
of  Markov-1  Toeplitz  matrices  has  been  shown  to  be  covered  by  the 
class  diagonal  in  the  Cosine  transform^-^. 

The  problem  of  evaluating  the  asymptotic  performance  of  a  trans¬ 
form  technique  is  formulated  as  follows:  For  a  given  transform  net 
does  the  diagonal  class  ^asymptotically  cover  y  T  y  ?  From  the 
definition  of  asymptotic  equivalence  of  nets,  and  the  invariance  of 
the  weak  norm  under  unitary  transformations  we  have: 
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0 


JD  3.  y  T  y*  iff  Y  Ei  T  • 


That  is,  our  problem  Is  to  determine  whether  the  signal  class  is 
asymptotically  covered  by  the  class  diagonal  in  y . 

Remark :  These  definitions  and  approach  is  a  slight  modifi¬ 
cation  of  the  one  presented  in  (Pearl  [10]). 

We  also  use  the  notation  A-*-*-  if  both  .A  covers  (3  and  (3 
covers  .A.  We  say  then  that  A  and  £3  are  asymptotically  equivalent. 

3.  DISCRETE  GENERATION  OF  TRANSFORMS  AND  THEIR  ASYMPTOTIC  PROPERTIES 


3 . 1  Numerical  Quadrature  and  Transform  Generation 

Numerical  quadrature  is  a  scheme  of  approximating  integrals  by 
finite  sums.  It  is  discussed  in  any  book  on  numerical  analysis  and 
to  a  deeper  extent  in  (Krylov  [11]).  Our  approach  is  to  study  the 
asymptotic  behavior  of  transform  nets  using  integral  approximations 
of  finite  sums.  A  brief  summary  of  numerical  quadrature  is  given 
i n  the  sequel . 

Let  X  be  a  real  interval  (usually  we  shall  take  X  =  [-1,1]  or 
[0,1]).  Let  dy(x)  be  a  measure  on  X.  Consider  integrals  of  the 
form  /  f(x)dp(x).  A  numerical  quadrature  is  a  sum  of  the  form 


qN(f> 


N-l  u  u 

l  f(«j ) 

i  =0 


N-l 


are  called  the  weights  of  the 


{x?}  the  discretization  points, 
1  i  =o 


N-l 


numerical  quadrature  and 
N  is  the  order  of  the  numerical  quadrature 


A  quadrature  is  exact 
of  order  N  if  for  any  polynomial  P(x)  of  order  <  N,  Q n ( P )  = 

/  P(x)du(x).  A  quadrature  converges  W.R.T.  a  given  class  of  func- 


X 

tions  F  if  for  any  f  CF  CL(f}—  /  f(x)dv(x). 
-  n  x 


Consider  the  case  where  dy(x)  is  a  probability  measure  over 


i 


x  =  [-1,1]. 


To  demonstrate  the  process  of  generating  unitary  trans 


forms  by  discretization,  we  present  the  Gauss-Jacobi  quadrature 

formula.  Let  (P^(x)}  be  the  orthonormal  polynomials  W.R.T.  dy(x) 

obtained  by  the  application  of  the  Grahm-Schmidt  process  to 

x^n  =  0,1, •••N.  It  is  known  (Szego  [12])  that  PN(x)  is  a  polynom- 

”  ] 

ial  of  order  N  which  possesses  N  distinct  roots  in  (-1,1);  {xN} 

1  i  -o 

The  Gauss-Jacobi  theorem  states  that  there  are  numbers  (Christoffel 

TFi  - — — x - 

Numbers)  {X.}  such  that 
i=o 

N  N-l  w 

(1)  Xj  >  0  for  i  =  0  ...  N-l,  l  xj  -  1. 

i=o  ' 


(2)  The  Gauss-Jacobi  quadrature  QN(f)  A  £  x^  f(x^),  is 
exact  of  order  2N-1 . 

(3)  Q^(f)  converges  for  any  f  for  which  /  f(x)dy(x)  converges. 


This  theorem  implies  the  following  procedure  for  generating 

orthogonal  transform  nets.  The  polynomial  P.  •  P.’  0  <  k,  j  <  w-1 

k  j 

is  of  degree  <  2N-1,  thus 


QN(Pk  *  Pj) 


,1 

Pk(x)Pj(x)dp(x) 

-1 


6 


kj’ 


therefore: 


N-l 

l 

i-o 


If  we  let  A  yly  Pj(x^)  then  is  a  unitary  transform 

net.  We  call  transform  generated  by  this  procedure  Gauss-Jacobi 
Transforms . 

Examples:  (1)  Let  X  =  [-1,1]  dy(x)  =  l==  then  Pn(x)  = 

ir/l-x2 

cos  mre  Vl,e'f,e  x  =  cosO)  are  the  corresponding  orthonormal  polynomials 


ISO 


The  G-J  transform  is  formed  by  discretizing  PQ  •••  Pp  ^  at  the 
zeros  of  Pp.  This  gives  the  cosine  transform  as  discussed  in 
(Ahmed  [8]). 


(2)  Let  X  =  [-it,,.]  du(x)  •  /TT'dx.P^x)  =  >e- 

where  x  =  cos0  are  the  Tchebychev  polynomials  of  the  second  kind 
and  give  rise  to  the  sine  transform 


<N  A 
sij  4 


s  i  n 


ITT 

n  +  1 


sir,  =  <*,  sin 


ITT 


where  are  normalization  factors. 


If  X  =  [-ir.ir]  and  if  we  consider  integrals  of  the  form 
1  i  0 

j  f(e  ')e  d0  then  we  have  similar  results  for  the  Newton- 


-ir 

Cotes  Quadrature  Scheme.  The  discretization  points  in  this  case 
are  equally  spaced  in  [0,  tt);  x^1  =  i  •  i  =  0  •  ••  N-l  and  {x^} 

are  selected  so  that  the  resulting  quadrature  is  exact  of  order 

4  p  A 

N-l.  If  we  take  the  system  (e  }p  of  orthonormal  functions  then 
the  exactness  of  the  quadrature  scheme  gives  the  familiar  Discrete 
Fourier  Transform. 

Using  arbitrary  measures  on  the  unit  circle  and  their  corres¬ 
ponding  orthonormal  trigonometric  polynomials  one  can  produce  many 
other  trigonometric  transforms. 

3 . 2  Termwise  Asymptotic  Behavior  Of  Transforms  -  Generalized 
Toeplitz  Matrices 

Let  X  =  [-1,  1],  dy(x)  a  probability  measure  on  X,  and  PjJx) 
the  corresponding  othonormal  polynomials.  If  Qp(f)  is  the  Gauss- 
Oacobi  quadrature  scheme,  then  by  G-J  theorem,  this  scheme  conver¬ 
ges  for  every  function  f,  such  that  /  f(x)  dy(x)  exists.  In  par- 
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f 


ticular,  if  f(x)  is  a  bounded  positive  function  then 


qN(pk f  V 

If  Dn  =  Diagonal  [f(xjj)  •  •  •  f(xfj  -j )]  and  CN  is  the  Gauss- Jacobi 
Transform  corresponding  to  { P ^ }  then  this  simply  implies  that 


[CN  Dn  cJ 


N  -*-® 


/  pk(x)Pj ( x ) f (x)dp(x) 


That  is,  if  we  look  at  the  diagonal  net  A  {DN}  and  the  G-J  trans 

form  net  y  A  {C^},  then  the  elements  of  the  matrices  of  the  trans- 

★ 

formed  net,  y  y  tend  asymptotically,  termwi se ,  to  the  elements 
of  the  fixed  (infinite)  matrix 

[M(f)]ki  =  /  Pk(x)Pi(x)f(x)dp(x). 

X 

This  means  that  the  class  diagonal  in  y  is  termwise  asymptotically 
equivalent  to  the  class  of  finite  sections  of  infinite  matrices  of 


This  motivates  a  further  examination  of  such  classes.  Given 
an  orthonormal  system  {$n(x)}n  on  the  measure  space  (X,  dy )  consi 
der  the  infinite  matrix  defined  by 


f(x)4>k(x)4»j.  (x)dy(x) 


where  f  is  positive  bounded  function.  With  each  such  matrix  we 
associate  a  net  of  finite  sections  of  M ( f )  (i,j  =  0,1,  -  *  *  N - 1 ) ;  we 
use  m^  to  denote  this  net.  The  class  of  all  such  nets  is  called 
Generalized  Toeplitz  class  associated  with 


The  class  diagonal  in  y  is  termwise  asymptotically 
equivalent  to  the  generalized  Toeplitz  class  asso¬ 
ciated  with  { P|,  (x  ) }  |  . 

Exampl e  1 :  The  generalized  Toepl i tz.  matrices  associated  with 
(^(x)  A  :  e^1T*cx  are  the  ordinary  Toeplitz  matrices  since: 

>TT 

,  CM(f)3kJ  -  ^  j  ei(k-J)x  f{x)dx  =  f(k-j) 

-ir 

(f(k-j)  is  the  k-j  Fourier  coefficient  of  f).  That  is,  M(f)  has 
constant  diagonals. 

From  here  on  we  designate  the  term  Toeplitz  matrix  to  matrices 
T  =  (t-jj)  of  the  form  t.j  *  t(|i-j|)  which  are  positive  definite. 
Such  matrices  are  the  autocovariance  matrices  of  stationary  signals. 

The  previous  result  implies  that  the  class  diagonal  in  the  DFT 
(circulant  matrices)  is  termwise  asymptotic  equivalent  to  the  Toep¬ 


litz  class.  This  parallels  the  result  of  Hamidi  and  Pearl 


[13] 


show¬ 


ing  the  vanishing  of  the  off-diagonal  elements  of  y  t  y* ,  with  y 
being  the  net  of  DFT  matrices. 

Example  2:  Consider  again  the  DCT  which  has  been  shown  to  be 

the  G-0  transform  associated  with  the  Chebychev  polynomials  of  the 

first  kind  P^x)  =  cos  nir0  0  =  cos~^x.  We  now  find  the  generalized 

Toeplitz  matrices  associated  with  this  orthogonal  system.  Let  f(x) 

be  a  positive  bounded  function  on  X  =  [-1,  1].  Then 
f 1  fir 


1 


[M<f)3kj  -  j 


i 


f  (e)cos  kite  cos  j  it  9  d8  =  27 


f (e )cos ( k-j )e  do 


+  27 


fTT 


f ( 0 ) cos ( k+ j )0  d0 
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M(f ) 


Thus,  a  typical  M(f)  is  the  sum  of  a  Toeplitz  matrix  and  a  Hankel 
matrix  (both  infinite),  with  the  same  parameters.  Loosely  stating, 
our  theorem  above  states  that  any  cosine-diagonal  net  tends  term- 
wise  to  such  a  generalized  cosine-Toepl itz  matrix.  If  we  can  prove 
that : 

(1)  The  generalized  cosine-Toepl i tz  class  is  an  asymptotic 
cover  of  the  Toeplitz  class 

and  (2)  Termwise  asymptotic  behavior  implies  ordinary  asymptotic 
equivalence  (in  weak  norm  sense), 

then  we  could  show  that  the  cosine  diagonal  class  covers  the  Toep¬ 
litz  class. 

In  summary,  starting  with  a  measure  space  (X  dp)  and  an  ortho¬ 
normal  system,  we  followed  two  paths  to  produce  two  matrix  classes: 

★ 

Y  y »  The  class  diagonal  in  a  corresponding  transform  net,  and 
m^,  a  generalized  Toeplitz  class.  Termwise,  these  classes  approxi¬ 
mate  each  other  asymptotically.  Faced  with  the  problem  of  showing 

* 

Y  A  Y  t  we  prefer  to  show 

(1 )  y  •&  Y  •* — *•'  my 

( 2 )  m^  * — '■+  t 

This  approach  is  motivated  by  the  fact  that  is  a  more  manageable 

* 

matrix  form  than  y  /S'  y  because  the  elements  of  any  net  in  m^  do 
not  vary  with  N.  Consequently  the  asymptotic  behavior  of  the  nets 
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in  are  more  transparent. 

3 . 3  From  Termwise  To  Weak  Equivalence 

Theorem  1  :  Let  X  =  [-1  1],  d ;.  ( x )  a  probability  measure  on  X 
which  has  a  density,  { P  ^  ( x  ) }  the  orthogonal  polynomials  correspond¬ 
ing  to  dp(x)  and  f  a  bounded  positive  function.  Define 


“ij  4  tCN  Dn  CN]ij 


(CN,  Dn  as  in  3.2) 


“ij  4  /  f(x)P1(x)Pj(xjd„(x) 

X 


N-  1 


Then  lim  ^  £ 

N  -►  <»  i  ,j=o 


N 

aij  "  aij 


=  0 


This  means  that  y  y  is  asymptotic  equivalence  to  m^. 

N 

Proof :  By  3.2  we  have  termwise  convergence,  i.e.  a,- - ->a-iV 

•  J  ^  ->cv>  *  J  ■  o 


Also,  from  [13]  (Theorem  6.2): 
N-l 


1  y 
N  i.3-0 


N-l 


i  ,J=o 


N-l 

N°w  l  l 


i  ,J  =  o 


N 

2 

ai  J 

| a  i  j 

2 

.  . 

N 

aij 

2 

aij 

-1 

1 

-1 


f  (x)dy(x) 


f^x)dy(x) 


1  N;]  r/  N  x2  .  2  1  2  N;]  N 

N  I  ^aii^  ai  j  ^  '  N  .  I  _aij  ai  j 


1  ,  J  =  o  J 


i  ,  j  =  o 


,1 


Let  N  -*■  so  then  the  first  sum  goes  to  2 


f ^ ( x ) dy (x ) . 


Now 


1  V  N  1  M  r  ,  s 


N-l 


1  ,j  =  o 


-1 

N-l 


i. 


N  /  N  x 

a  * ' ' a  i  j  a  i  j  ' 
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<  1  Y  a..  aN.  -  a.  . 
“  N  l,j-o  ,J  u  15 


and  (as  generalized  Fourier  coefficient  of  f)  lim  ai i 

i  -*•  c» 


Let  e  >0  and  let  N  be  such  that  a..  <  e  for  i,j  _>  N-l.  If  we 
choose  N  >  N  large  enough  such  that  •  -  ai .  <  e  for  i ,j  £  N-l 

v  I  '  J  1  J  | 

we  get  ’ 


Hq-1  N 


,  N 

1  r  0 

r  ,  l  “u  “u 


\°u 


,  V  N 

i  r  0 

r —  }  a  •  •  a  .  .  -  a  .  • 

No  i,j=N-l  lJ  'J  '> 


i  V  i  No 

-  e  TT  •  ^  ai i  +  fT  I  ,  aij  -  aii 
No  1J  o  i , J  =  N - ! 


r  0 

l  - 


These  two  sums  are  .bounded  because 


•,  IK  -  I 

N~  I  a--  tends  to  f(x)dy(x)  ([13]  (Theorem  6.2))  and 
i,j=o  J 

N 

a.1!  - >  0  as  i  ,  j  -*  =°.  Thus 

I  J 

t 

,  N-l  N  ,  N-l  2 

N  .  I  “jj  “ij  •  if  .  I  «ij  <  E  constant 

1  ,3=0  J  J  1  ,J=0  J 

which  proves  that 


2  N;]  N 

N  “ij  aU 


2  i  (x )dy (x  ) 


Thus 


1  N"1  N 

»  1.1-0 


N  *» 


Q.E.D. 


SUMMARY: 


We  have  shown  that  for  measures  with  density 

y*  Y  >nY 

This  means  that  to  find  out  whether  Y  *  jar  y  -=>  y  it  is  enough 

(a)  To  find  the  form  of  generalized  ToSplitz  matrices  corres¬ 
ponding  to  dy (x  ) , 

(b)  Compare  those  infinite  matrices  with  signal  class  O',  and 
prove  their  equivalence. 

3 . 4  Appl i cations 

1.  We  first  use  this  method  to  demonstrate  that  the  cosine- 
diagonal  class  is  an  asymptotic  cover  of  a  useful  Toeplitz  sub¬ 
class.  Consider  the  class  with  Y  being  the  CDT  net.  As  shown 
before,  a  typical  net  is  of  the  form 


In e  want  to  show  that  the  Toeplitz  part  is  asymptotic  equivalent  to 
the  Toeplitz  plus  Hankel  net.  This  means  that  the  Hankel  net  must 
tend  (in  weak  norm)  to  zero,  i.e.,  we  want  to  show 

lim  l  [ao  +  2al  +  3a2  +  *  '  '  N  -  aN-l]  =  0 

N  -*■  00 

Clearly,  if  we  restrict  our  Toeplitz  class  to  satisfy  the  additional 

00  2 

smoothness  condition:  l  m  aj;  <  «  then  the  limit  above  approaches 

m=o 
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■ax 


0  like  N 


i 


We  conclude,  therefore,  that  the  class  of  stationary  signals 

satisfying  £  m  ®  is  asymptotically  covered  by  the  cosine 

m  =  o 

transform.  This  is  a  substantial  generalization  of  Hamidi  and 
Pearl's^  result.  It  implies,  for  instance,  that  the  DCT  is  asymp¬ 
totically  optimal  for  all  finite-order  Markov  signals.  Moreover, 
the  proof  of  this  statement  is  surprisingly  simple  involving  none 
of  the  tedious  calculations  used  in  [9]  for  Markov-1  processes. 

2.  The  same  procedure  can  be  applied  to  polynomials  on  the 

unit  circle  (trigonometric  polynomials)  which  yield  the  DFT  by 

★ 

discretization.  In  this  case  the  equivalence  y  jy  y  ■*-*■  tr  has  a 
particular  meaning:  the  left  hand  side  are  the  circulant  matrices 
and  the  right  hand  side  is  just  the  ordinary  Toeplitz  class.  Thus 
we  readily  obtain  the  result  that  the  circulant  class  is  asymptotic 
equivalent  to  the  Toeplitz  class. 

CONCLUSIONS 

A  method  for  studying  the  asymptotic  behavior  of  spectral 

transformations  was  developed  using  numerical  quadrature  theory. 

Using  this  method,  the  asymptotic  optimality  of  common  unitary 

transforms  can  be  tested  conveniently.  A  practical  new  result 

obtained  by  this  method  states  that  the  discrete  cosine  transform 

is  asymptotically  equivalent  to  the  Ka rhunen-Lo6ve  transform  of 

”2 

stationary  signals  satisfying  £  m  <  »  (e.g.  finite  order 

m=o 

Markov  processes). 
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REAL  TIME  TELEVISION  IMAGE  BANDWIDTH  REDUCTION 
USING  CHARGE  TRANSFER  DEVICES* 

H.  J  Whitehouse 
R.  W  Means 
E  H  Wrench 


Abstract 


Recent  advances  in  analog  semiconductor  technology  have  made  possible  the  direct  sensing  and  processing  of  television  images.  By 
combining  a  charge  transfer  device  (CTD)  imager  and  a  CTD  transversal  filter,  real  time  image  sensing  and  encoding  have  been  achieved 
with  low  power  integrated  circuits  so  that  digital  transmission  and  bit  rate  reduction  are  made  possible  using  differential  pulse  code 
modulation  (DPCM).  Good  mean  square  error  performance  and  freedom  from  DPCM  artifacts  are  possible  in  a  hybrid  intraframe  image 
encoder  The  hybrid  transform  encoder  performs  a  discrete  cosine  transform  (DCT)  on  each  line  of  the  television  image  as  it  is  scanned. 
This  compacts  the  variance  into  low  frequency  coefficients  and  the  DPCM  encodes  the  corresponding  DCT  coefficients  between  successive 
lutes.  Computer  simulation  of  this  hybrid  coding  technique  has  shown  good  performance  on  256  x  256  pixel  images  at  0.5  bits/pixel 
and  channel  bit  error  rates  of  10*-  An  experimental  system  using  a  low  resolution  General  Electric  100  x  100  charge  injection  device 
camera  and  a  Texas  Instruments  bucket  brigade  transversal  filter  as  part  of  the  DCT  processor  has  been  constructed  and  provides  good 
low  resolution  image  quality  at  I  bit/pixe!  and  bit  error  rates  of  1 0*3  A  high  resolution  vidicon  compatible  system  is  also  being 
constructed. 


Introduction 

Unitary  transforms  for  image  encoding  have  been  used  for  intraframe  encoding.1 1 1  In  addition,  these  techniques  may  also  be 
applied  to  interframe  and  multispectral  encoding  However,  all  unitary  transformations  are  information  preserving  and  no  bandwidth 
reduction  results  from  the  application  of  the  transform  to  the  image.  Instead,  the  transforms  redistribute  the  variance  associated  with 
each  picture  element  (pixel)  so  that  subsequent  to  the  transform,  basis  restricting  operations  on  the  transform  coefficients  will  result  in 
bandwidth  reduction.  Upon  reconstruction  of  the  original  image  from  the  basis  restricted  transform  coefficients,  a  degraded  version  of  the 
original  image  can  be  obtained.  Unfortunately,  the  interrelationship  between  the  type  of  transform,  the  form  of  the  noninvertible  oper¬ 
ation.  and  the  type  of  degradation  in  the  reconstructed  image  is  very  complicated  and  subjective  The  universally  used  analytic  criterion 
of  the  mean-square-error  is.  at  present,  the  best  compromise  technique  for  transform  comparison. 

For  the  particular  operation  of  basis  restnciion  by  truncation,  a  particularly  simple  interpretation  of  the  bandwidth  reduction  can 
be  made  The  transforms  may  be  viewed  as  a  variance  redistributing  operation  that  approximately  decorrelates  the  transform  coefficients 
while  transforming  the  variance  associated  with  each  picture  element  into  the  low-order  coefficients  of  the  transform.  Under  the  assump¬ 
tion  that  each  set  of  picture  elements  can  be  considered  as  a  sample  function  from  a  wide  sense  stationary  random  process  with  correlauon 
function  riT',  there  exists  an  optimum  discrete  transformation,  the  Karhunen-Loeve  transformation,  which  totally  deconelates  the  trans¬ 
form  coefficients  and  maximally  compacts  the  variance  to  the  low-order  coefficients.  All  other  transformations  can  be  compared  in  their 
performance  by  comparing  their  transform  coefficient  decorrelation  and  variance  compaction  with  this  optimum  transformation. 

This  intuitive  interpretation  can  be  made  rigorous  through  the  use  of  the  rate-distortion  criterion.!—  3)  ]t  has  been  found  from 
experience  that  the  closer  the  eigenvectors  of  the  transformation  are  to  the  eigenvectors  of  the  optimum  Karhunen-Loeve  transformation 
the  greater  the  variance  is  compacted  and  the  more  the  coefficients  can  be  truncated  while  maintaining  a  fixed  rate  distortion  or  mean- 
squar  e-error. 


Transform  Encoding 


Karhunen-Loeve  Transformation 

If  a  continuous  time  function  of  zero  mean  and  autocorrelation  function  R(r)  *  e"alTl'  is  considered  to  be  a  sample  function  from  a 
wide-sense  stationary  random  process,  then  this  time  function  can  be  explicitly  expanded  by  the  Karhunen-Loeve  expansion!4*  and  the 
resulting  coefficients  will  be  uncorrelated.  For  a  discrete  function  of  zero  mean  and  autocorrelation  function  R(r)  *  t|r',  which  may  be 
considered  as  a  sample  function  from  a  first-order  Markov  process,  a  similar  discrete  Karhunen-Loeve  transformation  may  be  defined  *5* 
This  transformation  diagonalizes  the  covariance  matrix  and  is  optimal  in  the  mean-square-error  sense  for  a  restricted  set  of  basic  functions 
that  do  not  span  the  complete  space. 

The  discrete  Karhunen-Loeve  expansion  Is  given  by!5*  for  the  case  N  *  2m  as 

2m 

Gy  *  y.  - - — -  sin  ju>n  [k  -I2m  +  1 5/2)  +  n»/2  j  g„ 

n-I  -m  +  ^  (  ' 

k  *  1 . 2 . 2m  (I) 


where 


I  -r* 


I  -  2r  cos  wn  +  r* 


C) 


•To  be  published  in  5PIE  Conference  Proceedings  (. A  ug.  1975). 
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jnd  wn  are  the  first  N  positive  roots  of 


tan  2mu>  = 


-t  1  -r-)  sin  to 
(cos  to  -2r  +  r~  cos  to) 


(31 


Since  the  discrete  Karhunen-Loeve  expansion  involves  both  the  solution  of  a  transcendental  equation  and  the  evaluation  of  the  auto¬ 
correlation  function  of  the  data  to  be  transformed,  real  time  computation  of  this  transform  is  quite  complex  However.  Habibi  and 
Wintz1 1 )  have  shown  that  Karhunen-Loeve  transformations  calculated  using  approximate  autocorrelation  functions  are  satisfactory  for 
many  applications. 

Discrete  Fourier  Transform 

Since  the  discrete  Fourier  transform  is  asymptotic  to  the  Karhunen-Loeve  transformation^*  for  the  exponential  covariance  function 
and  the  basis  vectors  are  picture  independent,  the  Fourier  transform  represents  a  logical  choice  for  real  time  implementation.  The  Fourier 
transform  exists  for  all  data  lengths  N .  This  is  defined  by 


N-l 

Gk*£  e-itr2nk/N  ^  k  =  0.  1 _ N-l  (4) 

n=0 


Discrete  Cosine  Transform 


Two  different  types  of  discrete  cosine  transform  (DCT)  are  useful  for  reduced  redundancy  television  image  transmission.  Both  are 
obtained  by  extending  the  length  N  data  block  to  have  even  symmetry,  taking  the  discrete  Fourier  transform  iDFT)  of  the  extended  data 
block,  and  saving  N  terms  of  the  resulting  DFT  Since  the  DFT  of  a  real  even  sequence  is  a  real  even  sequence,  either  DCT  is  its  own  in¬ 
verse  if  a  normalized  DFT  is  used. 

The  "Odd  DCT”  (ODCT)  extends  the  length  N  data  block  to  length  2N  -  1 .  with  the  middle  point  of  the  extended  block  as  a  center 
of  even  symmetry  The  "Even  DCT”  (EDCT)  extends  the  length  N  data  block  to  length  3N.  with  a  center  of  even  symmetry  located  be¬ 
tween  the  two  points  nearest  the  middle  For  example,  the  odd  length  extension  of  the  sequence  A  B  C  is  C  B  A  B  C.  and  even  length  is 
C  B  A  A  B  C.  In  both  cases,  the  symmetrization  eliminates  the  jumps  in  the  periodic  extension  of  the  data  block  which  would  occur  if 
one  edge  of  the  data  block  had  a  high  value  and  the  other  edge  had  a  low  value,  in  effect,  it  performs  a  sort  of  smoothing  operation  with 
no  loss  of  information.  It  will  be  noted  that  the  terms  "odd"  and  "even”  in  ODCT  and  EDCT  refer  only  to  the  extended  data  block  - 
in  both  cases  the  extended  data  block  has  even  symmetry. 

Let  the  data  sequence  be  gg,  g| ,  .  .  .  ,  g>j_|  The  ODCT  of  g  is  defined  as 

N-l  -)2<rnk 

Gk=  gn  e  :N  *'  fork  =  0.  1 . N-l  (5) 

n=-i  N-l) 

where 


g-n  =  *n  f°m  =  0.  1 - N-l. 

By  straightforward  substitution  it  may  be  shown  that 

N-l  -i2trnk 
Gk  =  2  Re  Sj  e  ;N  *  1 
n=0 


where  2"is  defined  by  equation  (8). 


1 0.5  g^  n  =  0  ) 

\  1 
'  gp,  n  =  1 . N  -  1  * 

The  EDCT  of  g  is  defined  by  equation  (9),  where  the  extended  sequence  is  defined  by  equation  1 10). 

-ittk  n_|  -i2trnk 

"’N  X-  2N 

Gk  *  e  2.  *n  e  for  k  =  0.  1 . N-l 

n»-N 
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(9) 


*-l  -n"*n  for  n*  0,1 - N-l 


(10) 
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If  (he  mutually  complex  conjugate  terms  in  equation  <d)  are  combined,  then  equation  t 
an  alternate  way  of  defining  the  fcDCT. 


1 1 )  results  Equation  (III  may  be  viewed  as 
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Ahmed1 ' '  has  investigated  the  use  of  the  EDCT  as  a  substitute  lot  the  Kurhunen-Loexe  transform  and  linds  that  it  is  superior  to  the 
Fourier  transform  and  is  comparable  to  the  Karhunen-Loeve  i  K-L)  in  rate-distortion  pertormance  while  maintaining  the  computation 
simplicity  of  a  transform  which  does  not  depend  on  the  picture  statistics  Habibi1^ 1  has  shown  by  simulation  that  the  DCT  is  equivalent 
in  a  mean-square-erfor  sense  to  the  K-L  transform  under  basis  restriction 


Hybrid  Transforms 

Mixed  transforms  can  be  used  either  for  intraframe  encoding  or  mterframe  encoding  depending  on  the  available  memory  and  the 
type  of  transforms  implemented.  In  order  for  a  mixed  transform  to  be  competitive  with  one  of  the  conventional  two-dimensional  trans¬ 
forms.  it  must  offer  either  superior  performance  or  simplicity  of  implementation  Of  the  transforms  examined  the  odd  length  cosine 
transform  is  competitive  in  performance  since  it  can  b  plemented  as  the  real  part  ot  i  CZT  and  since  the  transform  samples  are  real 

and  are  the  samples  of  an  autocorrelation  function  wh  may  then  be  extrapolated  bx  well-known  techniques  In  simulations  of  trans¬ 

form  performance,  the  cosine  transform  has  been  shown  to  closely  approximate  the  behavior  of  the  Karhunen-Loeve  translorm 

The  benefits  of  mixed  transformation  implementation  and  minimum  memory  may  be  achieved  for  digital  transmission  by  combining 
a  one-  or  two-dimensional  unitary  transform  with  generalized  DPCM  in  a  hybrid  system  The  basis  of  operation  of  the  hybrid  transform 
is  that  the  unitary  transform  decorrelates  the  image  within  its  constraints  of  transform  type,  dimensionality,  and  block  size,  while  the 
generalized  DPCM  removes  the  correlation  between  transform  blocks.  This  hybrid  system  is  particularly  attractive  for  remote  sensor 
application  since  it  has  been  found  that  its  performance  is  approximately  as  good  as  the  Karhunen-Loeve  transform  and  its  implementation 
requires  minimum  memory  *8) 


Implementation 


Computational  Modules 

A  linear  transform  on  sampled  data  of  finite  extent  may  be  viewed  as  the  multiplicat  on  of  a  vector  by  a  matrix  Multiplication  by 
diagonal,  circulant,  or  Toeplitz  matrices  may  be  accomplished  rapidly  with  simple  computational  hardware  modules  Multiplication  by 
an  N  X  N  diagonal  matrix  requires  only  a  scalar  multiplier  and  a  memory  containing  N  values  to  provide  serial  access  to  the  reference 
function  Multiplication  by  an  N  X  N  Toeplitz  matrix  corresponds  to  a  convolution  and  may  be  performed  using  a  transversal  filter 
having  2N-I  taps  Multiplication  by  an  N  X  N  circulant  is  a  special  case  of  multiplication  by  N  X  N  Toeplitz  matrix  in  which  the  length 
of  the  transversal  filter  may  be  reduced  to  N  taps  if  the  data  block  is  recirculated  through  the  filter  or  reread  into  the  filter  from  a  buffer 
memory 

One-Dimensional  DFT 


Linear  filters  have  been  used  for  many  years  for  the  calculation  of  the  power  spectra  of  continuous  signals.  One  of  the  earliest 
methods  used  a  bank  of  wave  Filters  to  measure  the  spectra  in  fractional  octave  bands  for  telephone  network  equalization.^)  However, 
when  tncreased  resolution  was  required  the  number  of  filters  rapidly  become  unmanageable.  An  alternative  which  overcame  the  difficulty 
of  a  large  number  of  filters  each  with  small  time-bandwidth  product  was  to  substitute  one  linear  fm  (chirp)  filter  with  large  time-bandwidth 
product  and  to  employ  matched  filtering  In  this  system  the  signal  to  be  analyzed  is  used  to  single  sideband  (SSB)  modulate  a  locally 
generated  chirp  signal  and  the  composite  modulated  signal  is  filtered  in  a  chirp  delay  line  matched  filter.  Each  component  of  the  input 
signal  spectrum  shifts  the  locally  generated  chirp  to  a  different  position  in  the  spectrum  after  SSB  modulation  and  these  shifted  chirps 
then  correlate  with  the  reference  signal  represented  as  the  impulse  response  of  the  matched  filter  at  different  times.  Thus  the  output 
signal  amplitude-time  history  reflects  the  amplitude-frequency  composition  of  the  input  signal. 

Bleustein1 1  ®)  recognized  that  the  discrete  Fourier  transform  (DFT)  of  sampled  data  was  amenable  to  a  similar  interpretation.  In 
addition  to  just  calculating  the  magnitude  of  the  Fourier  transform,  linear  filters  could  calculate  the  phase  and  thus  all  of  the  operations 
such  a  cross  convolution  and  a  cross  correlation  could  be  calculated.  This  technique  came  to  be  called  the  chirp-Z  transform  tCZT)  and 
can  be  applied  to  other  problems  besides  just  the  calculation  of  the  DFT.’  *  *  *  Prior  to  these  developments,  digital  computation  of  the 
DFT  had  been  significantly  improved  by  the  use  of  a  special  algorithm  called  the  fast  Fourier  transform  (FFT)  which  was  described  by 
Cooley  and  Turkey  * '  -)  The  FFT  algorithm  gained  rapid  popularity  in  signal  processing  since  it  allowed  the  calculation  of  the  DFT  to 
be  done  using  significantly  fewer  machine  operations  (multiplications)  than  direct  evaluation. 

By  direct  inspection  it  is  observed  that,  if  symmetries  of  the  function  exp  j*2nm/N  are  not  exploited,  then  the  number  of  complex 
multiplications  required  will  be  N-  corresponding  to  N  multiplications  for  each  frequency  component  evaluated.  Even  on  high  speed 
digital  computers  this  can  become  the  limiting  consideration  in  signal  processing  applications.  The  advantage  of  the  FFT  algorithm  is 
that  for  highly  composite  values  of  the  DFT  size  N  the  number  of  multiplications  is  proportional  to  N  logsN. 

Although  the  FFT  has  been  successful  in  substantially  reducing  the  computing  time  and  cost  of  using  general  purpose  digital  com¬ 
puters  it  has  several  disadvantages  for  special  purpose  real  time  computation.  At  high  throughput  rates  which  art  required  for  real  time 
image  processing  the  processor  either  must  operate  logjN  times  faster  than  the  data  rate  or  pipeline  structures  which  use  distributed 
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memory  and  logiN  multipliers  must  be  used.  In  addition,  the  internal  arithmetic  of  the  FFT  processor  must  be  done  at  increased  preci¬ 
sion  in  order  to  compensate  for  the  multiple  round  off  errors  introduced  by  the  successive  stages  in  the  FFT  processor  Although  these 
difficulties  can  be  overcome,  it  is  not  always  possible  to  arrange  the  computation  in  a  form  where  the  sue  of  the  transform  is  highly 
composite.  For  the  above  reasons  and  because  of  the  difficulty  of  obtaining  small,  low  power,  fast  analog  to  digital  converters,  linear 
transversal  filter  implementations  of  the  chirp-Z-transform  are  attractive'  *  3>  rather  than  the  previous  <  ZT  implementation  which  used 
an  FFT  to  perform  the  required  convolution 

The  DFT  may  be  easily  reduced  to  the  form  suitable  for  linear  filtering  by  the  substitution 


'  n-  -  in  -  mi*  ■ 


1 1:> 


which  changes  a  product  of  variables  into  a  difference  so  that 
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This  form  is  seen  to  be  equivalent  to  factoring  the  Fourier  mal.ix  F  into  the  product  of  three  matrices 

F  »  DTD  1 14| 

where  D  is  a  diagonal  matrix  with  elements  dnn  »  exp  t-j»n*  N I  and  T  is  a  Toeplitz  matrix  with  elements  tnm  »  exp  (jinn  -  ml*  N). 

The  CZT  algorithm  is  easily  implemented  by  transversal  filter  techniques  In  this  case  the  DFT  is  computed  by  premultiplication 
by  a  discrete  chirp,  convolution  with  a  discrete  chirp,  and  postmultiplication  by  a  discrete  chirp  Figure  1  shows  this  configuration 
However,  it  must  be  remembered  that  both  the  multiplications  and  convolutions  are  complex  and  a  suitable  representation  of  the  com¬ 
plex  numbers  must  be  used  One  representation  is  by  real  and  imaginary  part  Figure  2  shows  the  DFT  organized  as  a  CZT  and  imple¬ 
mented  with  parallel  conputation  of  the  real  and  imaginary  parts.  In  Figure  2  the  input  signal  is  represented  asg  *jg|  and  the  out¬ 
put  signal  is  represented  as  G  -  Gb  +jG|,  where  it  is  understood  that  g  *gn  n  «  0.  ,.N  -  I  and  G  *  Gn  n»Q.  .N-l 

In  order  to  determine  the  specific  form  of  the  transversal  filters  it  is  necessary  to  know  the  specific  value  of  N  When  N  is  odd  the 
Toeplitzmatnx  T  may  be  represented  as  a  transversal  filter  wiih  2N  -  1  complex  taps  h_([q_Gto  h^j_i  where  hn  -  W~n“  -  n  -  -iN-'.  ■ 
to  N-l ).  and  W  =  exp  i-j2ir/N),  The  required  convolution  has  oeen  implemented  with  the  general  transversal  filter  shown  in  Figure  ■ 
When  N  is  even,  it  can  be  shown  tha  Tn  m  3  T^^n  m  whe-e  the  subscripts  are  reduced  mod  N.  Thus  T  is  a  circulant  mairix  and 
can  be  implemented  with  a  recirculating  transversal  filter  js  shown  in  Figure  4  where  the  number  of  complex  taps  is  N  and  :he  tap  w  eights 
are:  h„  »  Wn'  -  n-0 . N-l 

Charge  Coupled  Devices 


CCDs  are  sampled  data  analog  circuits  which  can  be  fabricated  by  Metal  Oxide  Semiconductor  l MOSt  technology  as  LSI  compo¬ 
nents/'4*  As  such  they  are  directly  compatible  with  other  MOS  circuits  Current  CCD  transversal  filters  have  operated  as  video  devices 
with  sample  rates  up  to  5  MHz.  CCDs  operate  by  the  manipulation  of  injected  minority  earners  in  potential  wells  under  MOS  capacitors 
and  thus  behave  as  capacitive  reactances  with  low  power  dissipation.  However,  since  the  potential  wells  which  contain  the  minonty  car¬ 
riers  also  attract  thermally  generated  minonty  earners,  there  is  a  maximum  storage  time  for  the  analog  signal  which  depends  on  the  dark 
currant  associated  with  the  temperature  of  the  silicon  Under  normal  conditions  at  room  temperature,  dark  aments  are  tens  of  nAmps, 
cm*  and  storage  times  of  hundreds  of  milliseconds  can  be  achieved. 

There  are  many  ways  in  which  unidirectional  charge  transfer  can  be  achieved  The  first  developed  was  a  three-phase  clocking  struc¬ 
ture  which  is  illustrated  in  the  transversal  filter  of  Figure  5.  The  three  electrode  CCD  structure  is  planar,  much  like  the  SAW  devices,  and 
the  direction  of  charge  propagation  is  determined  by  the  sequence  of  potentials  applied  to  the  three  electrodes.  Unfortunately,  if  the 
minority  carriers  are  allowed  to  collect  at  the  semiconductor-oxide  boundary,  poor  charge  transfer  efficiency  will  result  due  to  minonty 
carriers  getting  caught  in  trapping  sites.  This  means  that  the  CCD  will  behave  nonlinearly  unless  there  is  sufficient  propagating  charge 
present  to  fill  all  of  the  traps.  By  biasing  the  operating  condition  of  the  CCD  so  that  about  10%  of  the  dynamic  range  is  used  for  the 
injection  of  a  “fat  zero,"  the  traps  are  kept  continuously  filled  and  the  device  has  over  a  60  dB  dynamic  range.  In  practice,  a  video  signal 
representing  the  signal  to  be  processed  is  added  to  a  fixed  bias  somewhat  larger  than  one-half  of  the  peak-to-peak  value  of  the  signal. 
Since  the  effective  storage  time  of  the  device  is  long  relative  to  the  time  required  to  execute  a  convolution,  CCDs  can  be  considered  to  be 


*  Denotes  either  convolution  or  circular  convolution 


Figure  1 .  Chirp-Z  Transform  Implementation  of  the  DFT 


Figure  2.  DFT  Via  CZT  Algorithm  with  Parallel  Implementation 
of  Complex  Arithmetic 


Figure  5.  Schematic  of  the  Sampling,  Weighting,  and  Summing  Operation 


interruptible  signal  processors  and  as  such  are  more  compatible  with  the  executive  control  required  for  signal  processing.  A  64  point  CCD 
filter  with  discrete  cosine  transform  sine  and  cosine  chirps  is  shown  in  Figure  6.  This  chip  was  developed  by  Texas  Instruments  for  the 
Naval  Undersea  Center  for  image  processing.  In  addition  to  the  four  DCT  filters,  four  DFT  filters,  a  Hilbert  transform  and  other  experi¬ 
mental  signal  processing  functions  were  also  implemented. 

Current  research  in  CCDs  is  directed  toward  improving  the  charge  transfer  efficiency  and  removing  the  requirement  of  continuous 
“fat  zero”  charge  injection  by  ion  unplantation  techniques  which  keep  the  minority  earners  away  from  the  semiconductor  oxide  bound¬ 
ary.  Ion  implantation  is  also  being  used  to  provide  asymmetric  potential  wells  so  that  simpler  two-phase  clocking  can  be  employed.  C  ar- 
rently  available  CCDs  have  500  stages  with  0.0999  transfer  efficiency  and  devices  with  up  to  2000  stages  are  planned. 

Another  charge  transfer  device  similar  to  the  CCD  is  the  Bucket  Brigade  Device  (BBDl  This  is  a  sequence  of  MOS  transistors  coupled 
together  by  diffusion  enhanced  Miller  capacitance.  Although  these  devices  do  not  operate  at  frequencies  as  high  as  CCDs,  they  have 
better  low  frequency  performance  since  they  include  active  devices.  A  CZT  has  been  implemented  with  two  BBD  chips.  Two  200  tap 
filters  are  implemented  on  each  chip:  one  a  discrete  cosine  and  the  other  a  discrete  sine  filter  The  BBD  chip  is  shown  in  Figure  1 .  The 
complex  chirp  used  in  the  premultiplier  and  a  typical  input  and  output  are  shown  in  Figure  S.  The  input  is  an  offset  cosine  wave  and  the 
output  shows  a  D  C.  component  plus  a  response  at  the  cosine  wave  frequency.  These  filters  can  operate  at  100  kHz  and  have  tap  accura¬ 
cies  better  than  I  With  careful  control  of  geometry,  both  BBD  and  CCD  filters  with  tap  accuracies  approaching  0.  IT  should  be  possible. 
This  chip  was  also  developed  by  Texas  Instruments  for  the  Naval  Undersea  Center. 
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Figure  6.  64  Point  CCD  Filters 
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Figure  "  Bucket  Brigade  Chirp  Filters 
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Figured.  BBD  Performance 


Discrete  Cosine  Transform 


Let  the  data  sequence  be  gQ.  g|  Sn-|  The  ODCT  of  g  is  defined  in  ( 5 )  as 


N-l  -)2irnk 
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The  identity  1 1 5 )  may  be  used  to  obtain  the  CZT  form  of  the  ODCT  shown  in  equation  ( I  hi. 
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The  block  diagram  of  the  ODCT  is  shown  in  Figure  9.  Since  only  real  inputs  and  outputs  are  required,  a  simplified  implementation  is 
possible  and  is  shown  in  Figure  1 0. 

A  corresponding  implementation  may  be  found  for  the  EDCT.  The  EDCT  of  g  is  defined  by  equation  (9),  where  the  extended  se¬ 
quence  is  defined  by  equation  ( 10). 
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Figure  9  ODCT  Block  Diagram 
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Figure  10.  ODCT  Expanded  Block  Diagram 


If  the  mutually  complex  terms  above  arc  combined  jnd  the  identity  (15)  used  the  C ZT  lorm  of  the  EDCT  becomes 
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Simultaneous  Computation  of  the  DFT  and  (In-  l»CT 

The  close  relationship  between  the  DFT  and  the  DCT  permits  the  use  of  common  modules  to  simultaneously  compute  both  trans¬ 
forms.  This  may  be  accomplished  most  simply  when  an  1  OCT  is  computed  using  DFT  modules.  The  sum  m  the  EDCT  defining  equation 
( I T  may  be  interpreted  as  a  length  IN  DFT  ol  the  extension  of  the  function  g  by  N  zeros.  This  leads  to  the  configuration  shown  in 
Figure  1 1 .  Alternatively,  if  the  odd  and  even  frequencies  in  the  zero-filled  DFT  are  considered  separately,  they  may  be  computed  using 
length  N  DFT  modules  as  shown  in  Figure  1 1 

System  Descriptions 


Two  hybrid  DCT/DPCM  bandwidth  reduction  systems  have  been  selected  for  construction  and  evaluation.  A  block  diagram  of  the 
systems  is  shown  in  Figure  13.  The  first  uses  a  slow  scan  image  sensor  and  a  Bucket  Brigade  Device  (BBD)  transform  implementation. 

The  second  uses  an  ordinary  vidicon  sensor  and  a  Charged  Coupled  Device  (CCD)  transform  implementation. 

In  the  BBD  system  a  100  x  100  pixel  solid  state  sensor  is  used.  The  nominal  horizontal  line  scan  time  is  one  millisecond.  The 
nominal  frame  rate  is  10  frames  per  second  which  can  be  displayed  without  flicker  through  the  use  of  a  scan  convener.  The  1  millisecond 
line  scan  time  was  chosen  in  order  to  match  the  sensor  to  the  BBD  filter  which  operates  at  a  clock  rate  of  100  kHz  with  good  charge  trans¬ 
fer  efficiency.  At  10  frames  per  second,  image  motion  should  be  reproduced  well  enough  for  many  applications,  even  though  some  pic¬ 
ture  detail  is  lost  because  of  the  low  spatial  sampling  afforded  by  the  1 00  x  100  pixel  format.  Minimum  overall  bit  rate  is  achieved  by  a 
combination  of  zonal  filtering  and  variable  bit  assignment  with  low  spatial  frequencies  assigned  more  bits  of  quantization  than  high  spatial 
frequencies.' '  * 

The  second  bandwidth  reduction  system  is  compatible  with  a  standard  vidicon  camera.  It  uses  CCD  filters  for  the  cosine  transform 
which  will  operate  at  4.8  MH2  sampling  rate  with  a  block  size  of  32  pixels.  Compatibility  with  standard  television  fomw  ■»- --nmeH 
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Figure  1  I .  Computation  of  the  DFT  and  EDCT  Using  a  Single  Length  IN  DFT  Module 
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Figure  1 1.  Computation  of  the  DFT  and  EDCT  Using  Two  Length  N  DFT  Modules 
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Figure  13  Image  Transmission  System 
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in  as  many  aspects  as  possible  If  the  interlace  field  is  used  directly  as  the  input  to  the  transform  hard*  are.  a  resolution  of  approximately 
-40  lines  by  256  pixels  is  possible  at  fill  fields  sec.  This  is  equivalent  to  a  video  bandwidth  slightly  less  than  '  5  MHr  Computer  simula¬ 
tion  of  this  system  is  shown  m  Figure  14  lor  1  bit  pixel.  P  is  the  channel  bit  error  rate 

The  implementation  of  selected  transforms  via  transversal  filters  was  proposed  at  the  All  Applications  Digital  Computer  Conference.' '  3  > 
The  computation  of  the  discrete  Fourier  transform  has  been  demonstrated  using  surface  wave  devices.1  '  The  discrete  cosine  transtorm 
has  been  demonstrated  using  bucket  brigade  devices.' 

The  DOT  implementation  block  diagram  is  shown  in  Figure  10.  The  convolutions  are  performed  in  both  television  systems  by  charge 
transfer  devices  built  by  Texas  Instruments.  The  multiplications  are  performed  by  conventional  circuitry.  The  reference  functions  used 
in  the  pre-  and  post-multiply  can  be  stored  in  a  real  only  memory. 

Figure  15  is  the  implementation  block  diagram  for  the  DPCM  part  of  the  system.  The  memory  is  a  line  store  which  is  used  as  the  pre¬ 
dictive  element  in  the  DPCM.  The  cosine  transform  coefficients  stored  in  the  memory  are  subtracted  from  the  new  transform  coefficients 
as  they  enter  the  DPCM.  The  difference  is  quantized  and  transmitted  with  a  selectable  number  of  bits  per  coefficient,  the  number  being 
a  lunction  of  the  assumed  variance  of  the  coefficient.  The  difference  coefficients  are  then  added  to  the  previous  line  coefficients  to  create 
the  predictive  element  for  the  next  line. 

The  BBD  slow  scan  system  is  shown  in  Figure  16.  It  consists  of  a  discrete  cosine  transform,  a  DPCM.  a  channel  simulator  in  which 
bit  errors  can  be  injected,  an  inverse  DPCM,  and  an  inverse  DCT.  The  system  is  built  on  eight  wire  wrapped  boards.  A  blow-up  of  one  is 
shown  in  Figure  1 7.  The  two  chips  in  the  center  are  the  BBD  devices.  System  performance  is  illustrated  in  Figure  18  at  2  bits/pixel. 

Because  the  original  picture  has  only  100  x  100  pixels,  adjacent  pixels  are  not  as  well  correlated  as  in  the  256  x  256  simulation  and  band¬ 
width  compression  algorithms  do  not  work  as  well.  The  system  performance  at  one  bit  per  pixel  is  shown  in  Figure  19.  Results  from  the 
CCD  high  resolution  system  are  not  yet  available. 


Figure  14.  Computer  Simulation  of  CCD  System 
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Conclusion 


The  use  of  an  intraframe  hybrid  transformation  composed  of  a  horizontal  unitary  discrete  cosine  transform  and  a  vertical  first -order 
DPCM  has  been  shown  through  simulation  to  have  performance  closely  approximating  that  of  a  two-dimensional  Karhenun-Loeve  trans¬ 
form.  This  hybrid  transform  has  been  computed  in  real  time  with  minimum  complexity  and  memory  through  the  use  of  LSI  bucket 
brigade  or  charge  coupled  devices  and  conventional  digital  DPCM  implementation. 
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ABSTRACT 

Iranslorm  coding  techniques  have  been  explored  extensively  in  theoretical 
studies  and  by  simulation.  It  has  been  shown  that  a  significant  bandwidth  reduction 
can  be  achieved  in  many  applications  with  minimal  image  degradation  and  relative 
tolerance  to  channel  errors.  The  major  drawback  of  transform  image  coding  for 
real  time  television  applications  in  the  past  has  been  computational  complexity. 
Transforms  and  algorithms  have  been  demonstrated  w  ith  recently  developed  charge 
coupled  device  and  acoustic  surface  wave  technologies  which  makes  nearly  optimum 
image  transform  coding  feasible  at  real  time  television  rates. 


INTRODUCTION 

Unitary  transforms  for  image  encoding  have  been  used  for  intraframe  encoding.  * 
In  addition,  these  techniques  may  also  be  applied  to  interframe  and  multispectral 
encoding.  However,  all  unitary  transformations  are  information  preserving  and  no 
bandwidth  reduction  results  from  the  application  of  the  transform  to  the  image.  In¬ 
stead.  the  transforms  redistribute  the  variance  associated  with  each  picture  element 
(pixel)  so  that  subsequent  to  the  transform,  basis  restricting  operations  on  the  trans¬ 
form  coefficients  will  result  in  bandwidth  reduction.  Upon  reconstruction  of  the 
original  image  from  the  basis  restricted  transform  coefficients,  a  degraded  version  of 
the  original  image  can  be  obtained.  Unfortunately,  the  interrelationship  between 
the  type  of  transform,  the  form  of  the  noninvertible  operation,  and  the  type  of 
degradation  in  the  reconstructed  image  is  very  complicated  and  subjective.  The  uni¬ 
versally  used  analytic  criterion  of  the  mean-square-error  is.  at  present,  the  best  com¬ 
promise  technique  for  transform  comparison. 

l  or  the  particular  operation  of  basis  restriction  by  truncation,  a  particularly  sim¬ 
ple  interpretation  of  the  bandwidth  reduction  can  be  made.  The  transforms  may  be 
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viewed  as  a  variance  redistributing  operation  that  approximately  decorrelates  the 
transform  coefficients  while  transforming  the  variance  associated  with  each  picture 
element  into  the  low-order  coefficients  of  the  transform  Under  the  assumption  that 
each  set  of  picture  elements  can  be  considered  as  a  s  imple  function  from  a  w  ide 
sense  stationary  random  process  with  correlation  function  r  T|.  there  exists  an  opti¬ 
mum  discrete  transformation,  the  Karhunen-Loeve  transformation,  which  totally 
decorrelates  the  transform  coefficients  and  maximally  compacts  the  variance  to  the 
low-order  coefficients.  All  other  transformations  can  be  compared  in  their  perform¬ 
ance  by  comparing  their  transform  coefficient  decorrelation  and  variance  compaction 
with  this  optimum  transformation. 

This  intuitive  interpretation  can  be  made  rigorous  through  the  use  of  the  rate- 
distortion  criterion-.  It  has  been  found  from  experience  that  the  closer  the  eigen¬ 
vectors  of  the  transformation  are  to  the  eigenvectors  of  the  optimum  Karhunen- 
Loeve  transformation  the  greater  the  variance  is  compacted  and  the  more  the 
coefficients  can  be  truncated  while  maintaining  a  fixed  rate  distortion  or  mean- 
square-error. 

The  use  of  two-dimension  transforms  can  provide  improved  performance  over  the 
use  of  transformations  on  a  line-by-line  basis  '.  The  most  direct  approach  is  to  seek 
a  two-dimensional  transform  which  simultaneously  decorrelates  the  transform  co¬ 
efficients  and  compacts  the  variance  into  a  corner  of  the  two-dimensional  transform 
coefficient  space.  One  method  is  to  find  a  two-dimensional  transform  which  can  be 
represented  as  the  product  of  a  transform  in  one  direction  and  a  transform  in  the 
other  direction.  Assuming  that  a  two-dimensional  picture  can  be  considered  as  a 
sample  function  from  a  random  process  with  two-dimensional  correlation  r,  IT<  1  i:  !7-' 
i.e.,  with  a  correlation  coefficient  r,  in  direction  one  and  a  correlation  coefficient  r, 
in  direction  two,  then  the  optimum  discrete  transformation  is  the  successive  use  of 
two  Karhenun-L.oeve  transformations;  the  first  with  parameter  r,  .  and  the  second 
with  parameter  r2 

Also  of  interest  in  transform  encoding  is  block  size.  Lor  a  one-dimensional  signal 
the  block  size  is  the  number  of  elements  of  the  transform,  and  the  performance  of  the 
transform  improves  monotonically  with  increasing  block  size,  l  or  two-dimensional 
images,  transform  performance  also  increases  with  increased  number  ol  elements  in 
each  dimension  of  the  transform.  However,  two  dimensional  transforms  usually  re¬ 
quire  intermediate  memory  to  store  the  transform  coefficients  in  the  first  direction 
while  the  transform  is  being  computed  in  the  second  direction. 

Alternatively,  two-dimensional  transforms  can  be  mixed  transforms,  i.e..  differ¬ 
ent  horizontal  and  vertical  transforms.  Performance  increases  with  the  number  of 
elements  in  each  direction  of  the  transform.  However,  for  a  fixed  first  transform 
size,  memory  requirements  tend  to  increase  linearly  with  the  number  of  elements 
in  the  second  transform  direction  since  all  of  the  coefficients  must  be  stored  from 
the  first  transform.  The  amount  of  intermediate  memory  is  minimized  by  the  use 
of  a  small  block  size  for  the  image  in  the  second  direction,  but  performance  depends 
critically  on  the  choice  of  the  second  transform.  Thus,  the  choice  of  a  mixed  trans¬ 
form  interacts  with  the  overall  system  design  and  the  available  memory  for  coeffi¬ 
cient  storage. 
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Karhunen-Loeve  transformation 

If  a  continuous  tutu*  function  of  zero  mean  and  autocorrelation  function  R(r)  = 
c-olr|  ls  considered  to  lie  a  sample  function  from  a  vv  id  e-sense  stationary  random 
process,  then  this  time  function  can  he  explicitly  expanded  by  the  Karhunen-Loeve 
expansion-*  and  the  resulting  coefficients  will  he  uncorrelated.  For  a  discrete  func¬ 
tion  ol  zero  mean  and  autocorrelation  function  RitI  =  FT  .  which  may  he  considered 
as  a  sample  function  from  a  first-order  Markov  process,  a  similar  discrete  karhunen- 
1  oeve  transformation  may  he  defined.-'  T  his  transformation  diagonalizes  the  co- 
variance  matrix  and  is  optimal  in  the  mean-square-error  sense  tor  a  restricted  set  of 
basic  functions  that  do  not  span  the  complete  space. 

flic  discrete  karhunen-1  oeve  expansion  is  given  by-'  for  the  case  N  =  2m  as 


2m 


( ii.  - 


2  - - - -  sm  /con  |k  -  (2m  +  I )  2|  +  nrr  :t  gn 
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where 


V  = 


1  -r- 


I  -  2r  cos  +  r- 
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and  can  are  the  first  N  positive  roots  of 


tan  2moj 


1 

-(I  -  r~)  sin  u 
(cos  u;  -  2r  +  t~  cos  co) 


(3) 


Since  the  discrete  Karhunen-Loeve  expansion  involves  both  the  solution  of  a  trans¬ 
cendental  equation  and  the  evaluation  of  the  autocorrelation  function  of  the  data  to 
be  transformed,  real  time  computation  of  this  transform  is  quite  complex.  However. 
Habibi  and  Wintz*  have  shown  that  Karhunen-Loeve  transformations  calculated  using 
approximate  autocorrelation  functions  are  satisfactory  for  many  applications. 
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Discrete  Fourier  transform 


Since  the  discrete  i  ourier  transform  is  asymptotic  to  the  k  itluiiu  n-l  oevo  tint' 
formation1'  for  the  exponential  covariance  function  and  the  hav'  'cciois  arc  picluic 
independent.  the  Fourier  transform  represents  a  logical  choice  lor  icul  time  iinplcmcn 
tation.  Hie  Fourier  transform  exists  for  all  data  lengths  N  It  isdclincu  In 


N-l 

r,k  =  2  e"i7r:nk/N  gM  k  =0.  1.  \  -  I 

n=() 


Many  methods  exist  for  the  computation  of  discrete  Fourier  coefficients.  I  lie 
( loertzel  algorithm  requires  a  number  of  computations  proportional  to  V  lu.i  can 
he  used  for  all  lengths  N.  When  N  is  highly  composite  "last  translormatioits  <.an 
be  used.7  Thus,  if  N  is  of  the  form  21'.  then  the  number  ol  computations  can  he 
made  proportional  to  Nq.  Although  "last"  algorithms  have  been  successtullv  list'd 
tin  general  purpose  computers,  they  are  too  slow  for  real  tune  compulation  since 
the  algorithm  iterates  q  times  bet  ore  achieving  a  solution  I  his  problem  can  be  over¬ 
come  by  the  use  of  q  processors  in  a  pipeline  architecture.*  although  tins  increases 
the  complexity  of  the  processor. 

A  linear  filter  implementation  also  exists  for  the  discrete  Fourier  transform  which 
is  both  easily  implemented  and  suitable  for  real-time  computation  I  Im  algonthm 
vailed  the  clurp-Z  transform*1  is  based  on  the  substitution  2nk  -  n:  +  h  -  in  -  kr 
and  can  be  used  for  any  length  sequence  N.  With  this  substitution  the  Dl  I  becomes 

N-l 

=  e-ijrk~/N  2  citr(n-k)-/N  e-i7rn-/N  ,C| 

n=0 


The  transform  may  be  performed  as  a  premultiplicalion  by  a  discrete  chirp,  convolu¬ 
tion  with  a  discrete  chirp  of  twice  the  length,  and  postmultiplicalion  by  a  disciete 
chirp.  This  convolution  may  be  computed  with  either  acoustic  surface  wave  lilters  or 
charge  transfer  devices. '  0 

Discrete  cosine  transform 

Two  different  types  of  discrete  cosine  transform  (DCT)  are  useful  for  reduced  re¬ 
dundancy  television  image  transmission.  Both  are  obtained  by  extending  the  length 
N  data  block  to  have  even  symmetry,  taking  the  discrete  Fourier  transform  <  DT  F )  ol 
the  extended  data  block,  and  saving  N  terms  of  the  resulting  DFT.  Since  the  DFT 
of  a  real  even  sequence  is  a  real  even  sequence,  either  DCT  is  its  own  inverse  if  a 
normalized  DFT  is  used. 

The  "Odd  DCT”  (ODCT)  extends  the  length  N  data  block  to  length  2N-1 .  w  ith 
the  middle  point  of  the  extended  block  as  a  center  of  even  symmetty.  The  "Even 


I  Ml"  il  IM  I  (extends  the  length  N  data  Mock  lo  length  2V  with  a  center  of  even 
symmetry  located  between  the  two  points  nearest  the  middle  For  example,  the  odd 
length  extension  ol  the  sequence  A  B  C  is  (l)AH  C.  and  the  even  length  is  C  B  A 
A  It  (  In  both  cases,  the  syminelri/ation  eliminates  the  jumps  in  the  periodic  ex¬ 
tension  of  the  data  block  which  would  occur  if  one  edge  of  the  data  block  had  a 
high  value  and  the  other  edge  had  a  low  value:  in  effect  it  performs  a  sort  of  smooth¬ 
ing  operation  with  no  loss  ol  information.  It  will  be  noted  that  the  terms  "odd" 
and  "even"  in  Ol X "I  and  I  IX  1  refer  only  to  the  length  of  the  extended  data  block 
in  both  cases  the  extended  data  block  has  even  symmetry. 

Both  types  ol  |)(  I  may  be  implemented  using  compact,  high  speed,  serial  access 
hardware,  in  structures  similar  to  those  previously  described  for  the  Chirp-/  trans- 
lorm  t(  /  I  i  implementation  ol  the  1)1  I 


let  the  data  sequence  be  g^.  g. .  .  g\_|  The  OIX1  of  g  is  defined  as 
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w  here 


K_n  '  b"  »  =  <)•  I-  •  *  "  I 
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By  siraighllorwanl  substitution  it  may  be  shown  that 
\_|  -i2irnk 


n=() 


where  g  is  defined  by  equation  (X). 
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fhe  identity  ( 10)  may  be  used  to  obtain  the  CZT  form  of  the  OIXT  shown  in 
equation  till 


2nk  =n*  ♦  k*  -  in  -  ki* 


(10) 
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The  E LXT  ot  g  is  defined  by  equation  (12).  where  the  extended  sequence  i-.  di 
fined  by  equation  (13) 


-iffk  N-l  -i2ffnk 


r  =  .  :n 

VII.  c 
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for  k  =  0.1 .  .  N  - 


( I ;» 


n=-N 


g_|  _n  =  gn  ,or  n  =  1 . N  ~  I 


(13) 


If  the  mutually  complex  conjugate  terms  in  equation  (12)  are  combined,  then 
equation  ( 14)  results.  Equation  ( 14)  may  be  viewed  as  an  alternate  way  of  defining 
the  HIX'T 


Equation  (14)  may  be  put  in  the  CZT  format  given  in  equation  (15) 
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Ahmed  has  investigated  the  use  of  the  EDC'T  as  a  substitute  for  the  Karlnmen- 
Loeve  transform  and  finds  that  it  is  superior  to  the  Fourier  transform  and  is  com¬ 
parable  to  the  Karhunen-Loeve  (K-L)  in  rate-distortion  performance  while  maintain¬ 
ing  the  computation  simplicity  of  a  transform  which  does  not  depend  on  the  picture 
statistics.  Habibi*  -  has  shown  by  simulation  that  the  DCT  is  equivalent  in  a  mean- 
square-error  sense  to  the  K-L  transform  under  basis  restriction. 
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I  igure  I  shows  the  performance  of  tile  two-dimciiNion.il  IH  I  on  .1  25<>x  2s<> 
sampled  picture  with  X  bit  quantization  P  is  the  hit  error  probability  of  the  trans- 
mission  channel. 


Hybrid  transforms 

Mixed  transforms  can  be  used  either  for  intraframe  encoding  or  interframe  encod¬ 
ing  depending  on  the  available  memory  and  the  type  of  transforms  implemented.  In 
order  for  a  mixed  transform  to  be  competitive  with  one  of  the  conventional  two- 
dimensional  transforms,  it  must  offer  either  superior  performance  or  simplicity  of 
implementation.  Of  the  transforms  examined  the  odd  length  cosine  transform  is 
competitive  in  performance  since  it  can  be  implemented  as  the  real  part  of  a  CZT 
and  since  the  transform  samples  are  real  and  are  the  samples  of  an  autocorrelation 
function  which  may  then  be  extrapolated  by  well-known  techniques.  In  simulations 
ot  transform  performance,  the  cosine  transform  has  been  shown  to  closely  approxi¬ 
mate  the  behavior  of  the  Karlnmcn-Foeve  transform. 

I'he  benefits  of  mixed  transformation  implementation  and  minimum  memory 
may  be  achieved  for  digital  transmission  by  combining  a  one-  or  two-dimensional  uni¬ 
tary  transform  with  generalized  DPCM  in  a  hybrid  system.  I'he  basis  of  operation  of 
the  hybrid  transform  is  that  the  unitary  transform  decorrelates  the  image  within  its 
constraints  of  transform  type,  dimensionality,  and  block  size,  while  the  generalized 
DPCM  removes  the  correlation  between  transform  blocks.  This  hybrid  system  is  par¬ 
ticularly  attractive  for  remote  sensor  application  since  it  has  been  found  that  its  per¬ 
formance  is  approximately  as  good  as  the  Karhunen-I.oeve  transform  and  its  imple¬ 
mentation  requires  minimum  memory.'  - 

Figure  2  shows  the  performance  of  a  hybrid  DC T/DPCM  simulation.  For  a  bit 
error  probability  of  P  =  0,  there  is  no  significant  difference  between  the  two- 
dimensional  DCT  and  the  DCT/ DPCM.  When  bit  errors  do  exist  in  the  transmission 
channel,  however,  the  hybrid  system  is  capable  of  better  performance.  This  is  shown 
for  a  bit  error  probability  P  =  It)--.  Figures  3.  4  and  5  show  the  simulation  of  the 
hybrid  DCT/DPCM  performance  lor  pictures  of  other  scenes  with  different  statistics. 


SYSTLM  IMPLEMENTATION 

Two  hybrid  DCT/DPCM  bandwidth  reduction  systems  have  been  selected  for 
construction  and  evaluation.  A  block  diagram  of  the  systems  is  shown  in  figure  6. 
The  first  uses  a  Charge  Injection  Device  (CID)'  image  sensor  and  a  Bucket  Brigade 
Device  (BBI)I  transform  implementation.  The  second  uses  an  ordinary  vidicon 
sensor  and  a  Charge  (('•pled  Device  (CCD)  transform  implementation. 

In  the  (  ID  system  a  100x100  pixel  solid  state  sensor  will  be  used.  The  nominal 
horizontal  line  scan  will  be  one  millisecond  The  nominal  frame  rate  will  be  10 
frames  per  second  which  can  be  displayed  without  flicker  through  the  use  of  a  scan 
converter.  The  I  millisecond  line  scan  time  was  chosen  in  order  to  match  the  sensor 
to  the  ttBD  filter  which  operates  at  a  clock  rate  of  100  kHz  with  good  charge  trans¬ 
fer  efficiency.  At  10  frames  per  second,  image  motion  should  be  reproduced  well 
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Cos/DPCM  0.5  bits/pixel,  P  =  0  Cos/DPCM  0.5  bits/pixel,  P  =  1 0-3  Cos/DPC 


Transmitter 


SOURCE  ENCODER 


DISPLAY 


Figure  6.  Image  Transmission  System 


enough  lor  many  applications,  even  though  some  picture  detail  will  he  lost  because 
ol  the  low  spatial  sampling  afforded  by  the  1 0U\1 01)  pixel  lorniat. 


Minimum  overall  hit  rate  will  be  achieved  In  a  combination  of  zonal  filtering  and 
variable  bit  assignment  with  low  spatial  frequencies  assigned  more  bits  of  quantization 
than  high  spatial  frequencies.  *  Table  I  shows  the  bit  rate  which  results  from  three 
overall  hits/pixel  assignments  at  a  pixel  rate  of  I  (P  pixels/see  An  overall  bit  assign¬ 
ment  ol  I  hit/pixel  should  result  in  a  signal-to-distortion  ratio  of  approximately  30 
dB.  In  addition,  since  channel  errors  occur  in  the  Fourier  domain,  channel  error 
rates  as  large  as  P  =  1 0_-  will  still  provide  useful  reconstructed  images. 


The  second  bandwidth  reduction  system  will  be  compatible  with  a  standard  vidi- 
con  camera.  It  will  use  (’CD  filters  for  the  cosine  transform  which  will  operate  at  a 
4.8  Mil/  sampling  rate  with  a  block  size  of  32  pixels.  Compatibility  with  standard 
television  format  will  he  maintained  in  as  many  aspects  as  possible.  If  the  interlace 
field  is  used  directly  as  the  input  to  the  transform  hardware,  a  resolution  of  approxi¬ 
mately  480  lines  by  256  pixels  is  possible  at  60  fields/sec.  T  his  is  equivalent  to  a 
video  bandwidth  slightly  less  than  2.5  Mil/.  Table  2  shows  the  bit  rate  which  re¬ 
sults  for  the  video  portion  of  the  field  . 


Table  I.  Bit  Rate  as  a  Function  of  (Quantization  for  the 
(  ID  System 
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I\i Me  2  Hit  Kate  as  a  I  unction  of  Quantization 
lor  the  Vidicom  Compatible  System 


Bits/Pixel 

Bit  Kate 
Megabits/sec 

*> 

0.6 

1 

4.X 

i/: 

2.4 

I  lie  implementation  of  selected  transforms  via  transversal  filters  was  proposed  at 
lire  All  Applications  P  gital  Computer  Conference.  *  ^  7  lie  computation  ol  the  dis¬ 
crete  I -ourier  transform  has  been  demonstrated  using  surface  wave  devices.  ^  The 
discrete  cosine  transform  has  been  demonstrated  using  bucket  brigade  devices.  ^ 

I  he  DC' I  implementation  block  diagram  is  shown  in  figure  7.  The  convolutions 
are  performed  in  both  television  systems  by  charge  transfer  devices  built  by  Texas 
Instruments.  The  multiplications  are  performed  by  conventional  circuitry.  Figure  X 
is  an  expanded  version  of  figure  7  showing  individual  components.  The  multipliers 
and  adders  can  be  put  on  the  same  chip  as  the  transversal  filters  so  that  the  complete 
cosine  transform  can  be  computed  using  two  chips.  The  reference  functions  used  m 
the  pre-  and  post-multiply  can  be  stored  in  a  read  only  memory. 

I  igure  '>  is  the  implementation  Mock  diagram  for  the  L)PCM  part  of  the  system. 

I  he  memory  is  a  line  store  which  is  used  as  the  predictive  element  in  the  DPCM. 
the  cosine  transform  coefficients  stored  in  the  memory  are  subtracted  from  the 
new  transform  coefficients  as  they  enter  the  DPCM.  The  difference  is  quantized  and 
transmitted  with  a  selectable  number  of  bits  per  coefficient,  the  number  being  a  func¬ 
tion  of  the  assumed  variance  of  the  coefficient.  The  difference  coefficients  are  then 
added  to  the  previous  line  coefficients  to  create  the  predictive  element  for  the  next 
line. 


-(INI  -  1)  <  m  <  (N  -1) 


0<n  <  (N  -  1) 


0<m<  (N  -  1) 


Figure  7.  Serial  Access  Implementation  of  the  DCT. 
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CONCLUSION 


The  use  of  an  mtratrame  hybrid  transformation  composed  of  a  horizontal  uni¬ 
tary  discrete  cosine  transform  and  a  vertical  first-order  DPCM  has  been  shown 
through  simulation  to  have  performance  closely  approximating  that  of  a  two- 
dimensional  Karhenun-Loeve  transform.  This  hybrid  transform  may  be  computed 
in  real  time  with  minimum  complexity  and  memory  through  the  use  of  LSI  bucket 
brigade  or  charge  coupled  devices  and  conventional  digital  DPCM  implementation. 


REFERENCES 

1 .  Habibi,  A.  and  Wintz,  P.A.,  “Image  coding  by  linear  transforma’ion  and  block 
quantization,”  IEEE  Trans  Commun.  Technol .  vol.  COM-19.  pp.  50-62.  Eeb. 

1971, 

2.  Tasto.  M.  and  Wintz,  P.A.,  “A  bound  on  the  rate-distortion  function  and  appli¬ 
cation  to  images,”  IEEE  Trans.  Inform.  Theory,  vol.  1T-I8.  pp.  150-150.  Jan. 

1072. 

3  Sakrison,  D.J..  and  Algazi,  V.R..  "Comparison  of  line-by-line  and  two-dimensional 
encoding  of  random  images,"  IEEE  Trans  Inform.  Theory,  vol.  IT-1  7.  pp.  386-308. 
July  1071. 

4.  Davenport,  W.B.,  and  Root,  W.L..  Random  Signals  and  Xoise,  McGraw-Hill: 

New  York  1058,  pp.  06-101 . 

5.  Ray,  W.D.,  and  Driver.  R.M..  “Further  decomposition  of  the  Karhunen-Loeve 
series  representation  of  a  stationary  random  process,"  IEEE  Trans.  Inform. 

Theory,  vol.  1T-I6,  pp.  663-668,  Nov.  1970. 

6.  Pearl,  J .,  "On  coding  and  filtering  stationary  signals  by  discrete  Fourier  trans¬ 
forms,”  IEEE  Trans.  Inform.  Theory,  vol.  IT-19,  pp.  220-232.  March  1073. 

7.  Cochran,  W.T.,  et  al.,  "What  is  the  Fast  Fourier  Transform,”  Proc.  IEEE,  vol. 

55,  pp.  1664-1674,  Oct.  1967. 

8.  Groginsky,  H.L.,  and  Works,  G.A.,  “A  pipeline  Fast  Fourier  Transform,”  IEEE 
Trans.  Computers,  vol.  C-19,  pp.  1015-1019,  Nov.  1970. 

9.  Rabiner,  L.R.,  et  al.,  “The  Chirp-Z  Transform  algorithm,”  IEEE  Trans  Audio 
and  Electroacoustics,  vol.  AU-17,  pp.  86-92,  June  1969. 

10.  Whitehouse,  H.J.,  Speiser,  J.M.,  and  Means,  R.W.,  "High  speed  serial  access  lin¬ 
ear  transform  implementations,”  presented  at  the  1073  Symp.  All  Applications 
Digital  Computer  (AADO,  Orlando,  Florida. 

1 1 .  Ahmed,  N.,  Natarajan,  T.,  and  Rao,  K.R.,  “On  image  processing  and  a  discrete 
cosine  transform,"  IEEE  Trans  on  Computers.  Vol.  C-23.  pp.  90-93,  Jan.  1974. 


193 


I  2  llabibi.  A..  “Hybrid  Coding  of  Pictorial  Data.  '  ( omm  f'ci  h  .  \o|  t  OM-22.  pp 
<»I4-o24.  Ma\ 

I  3  Michon.  C  J.  and  Burke.  H  V,  “Operational  characteristics  ol  l  II)  unagci. 
ISSCC Digest  of  Technical  Paper*,  pp.  2o-27.  Feb.  ll,74. 

14.  Alstip.  J.M..  Means.  R.W..  and  Whiteliouse.  H.J  ..  Real-Time  Discrete  homier 
Transforms  Using  Surface  Acoustic  H'<m’  Devices.  Proceedings  of  the  II  I  I  In¬ 
ternational  Specialist  Seminar  on  C  omponent  Performance  and  System  Appli¬ 
cations  of  Surface  Acoustic  Wave  Devices,  held  at  Aviemore.  Scotland.  24-2S 
Sept.  1 073. 

15.  Means,  R.W.,  Buss.  D.D.,  and  Wlnteliouse.  H.J..  Real-Time  Discrete  Fourier 
Transforms  Using  Charge  Transfer  Devices.  Proceedings  of  the  CCI)  Applications 
Conference  field  at  the  Naval  Hectronics  Laboratory  Center.  San  Diego.  Calif 
18-20  Sept.  1073.  pp.  05-101. 


194 


«  4H' 


I 


APPENDIX  K 

FREQUENCY  SYNTHESIS  V  A  THE  DISCRETE 
CHIRP  AND  PRIME  SEQUENCE  ROMS 


I  IvI.QlIliNCY  SYNTHESIS  VIA  THE  DISCRETE 
CHIRP  AND  PRIME  SEQUENCE  RONE 


James  M.  Alsup 
Harper  J.  Whitehouse 
Naval  Undersea  Center 
San  Diego*  CA  92132 


ABSTRACT:  Surface  acoustic  wave  sampled  data  filters  can  be  utilized 
as  serial  access  re ad -only -memories  to  directly  implement  at  carrier  freq¬ 
uencies  a  coherent  fast-frequency-hop  synthesizer  in  the  VHF  and  UHF  ranges. 
An  example  of  frequency  synthesis  using  surface  acoustic  wave  discrete  chirp 
filters  is  shown.  A  frequency  synthesis  scheme  using  a  surface  acoustic 
wave  prime- sequence  filter  is  described. 


To  be  published  in  the  Proceedings  of  the  IEEE*  Special  Issue  on 
Surface  Acoustic  Wave  Devices,  May  1976 
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FREQUENCY  SYNTHESIS  VIA  THE  DISCRETE 
CHIRP  AND  PRIME  SEQUENCE  ROMS 


Introduction. 

Many  types  of  frequency  synthesizers  have  been  developed  over  the  years  [1], 
but  only  recently  have  read-only-memory  (ROM)  synthesizers  begun  to  reach  maturity. 
This  has  occurred  primarily  because  of  the  advent  of  sine-cosine  lookup  tables  im¬ 
plemented  in  digital  hardware.  However,  other  types  of  ROM  devices  are  becoming 
available,  and  in  particular  surface  acoustic  wave  (SAW)  devices  are  starting  to 
serve  this  purpose.  ROM  synthesis  is  particularly  suited  for  applications  which 
require  coherent  fast  frequency  hopping.  SAW  ROM  technology  makes  possible  the 
direct  extension  of  ROM  synthesis  into  the  VHF  and  UHF  ranges. 

ROM  Synthesis. 

There  are  several  ways  to  synthesisze  a  sampled  sinusoid  using  discrete  ROMs. 

One  method  utilizes  one  period  of  a  sampled  sinusoid  stored  in  sequential  order 
in  a  random  access  ROM  [2,3].  If  all  samples  of  the  ROM  are  read  in  sequence  with 
a  sample  time  interval  At,  then  a  sinusoid  of  period  T  and  frequency  f0  is  generated, 
where  T  =  N  At  =  l/f0.  However,  if  every  kth  sample  is  read  with  the  same  sample 
interval  At,  then  the  frequency  of  the  generated  sinusoid  is  kf0.  Since  k  can  take 

on  the  values  0,1,..., N-l,  a  total  of  N  frequencies  which  are  all  harmonically  re¬ 
lated  can  be  generated.  The  sampling  ambiguity  called  aliasing  may  be  evident 

depending  upon  whether  complex  or  real  sinusoids  are  generated. 

A  second  method  for  the  generation  of  sampled  sinusoids  has  special  application 
in  SAW  technology  because  it  utilizes  serial  access  ROMs.  This  method  consists  in 
multiplying  the  outputs  of  two  N-sample  discrete  chirp  ROMs  whose  serial  readouts 
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are  mutually  processed  by  m  samples  to  synthesize  the  frequency  mf0.  For  complex 
ROMs  of  the  form  exp ( j  n  n^/N) ,  the  product  takes  the  form 

exp(j7rn^/N)  exp(-j  rr(n-m)^/N)  =  exp( j2 rrmn/N)  exp(-j  Trm2/N)  (1) 

where  n  corresponds  to  the  time  index  and  m  corresponds  to  the  frequency  index. 

The  right-hand  side  of  this  equation  can  be  interpreted  as  a  sampled  complex 
sinusoid  of  frequency  index  m  multiplied  by  a  complex  phase  shift  which  is  dependent 
only  on  m.  a1  so  note  that  expCjrrn^/N)  is  a  periodic  function  with  period  Nz)t. 

For  frequency  hopping  applications,  either  of  these  two  methods  requires  a 
hopping  interval  equal  to  or  greater  than  T  =  N  t  so  that  least  one  period  of 
each  frequency  is  synthesized  with  a  common  noivual  bandwidth  less  than  or  equal 
to  1/T.  Thus,  the  maximum  hopping  rate  for  these  typ>s  of  synthesizers  is  equal 

to  fQ.  Synthesis  may  take  place  at  a  carrier  frequency  fc,  m  which  case  f0  refers 
to  the  ''fundamental''  frequency  which  would  be  observed  if  the  comb-like  band  of  har¬ 
monics  t.o  be  generated  at  the  carrier  fc  were  shifted  to  baseband. 

synthesis  via  read-only-memory  is  attractive  because  it  makes  possible  the 
coherent  generation  of  a  harmonic,  group  of  discrete  frequencies.  These  frequencies 
all  will  start  with  known  initio!  phases,  and  us*.  uris  a  priori  information  ran  U 
made  to  build  coherent  receivers  capable  of  recognizing  particular  phase  relationshii 
among  sequentially-produced  tone  bursts  within  the  operating  band  of  the  synthesizer. 

ROM  synthesis  is  also  attractive  because  nearly  instantaneous  changes  can  be 
made  when  hopping  from  one  frequency  to  the  next.  In  contrast,  a  phase  lock  loop 
synthesizer  invariably  requires  several  periods  of  the  waveform  being  synthesized 
before  a  stable  lock  is  acquired. 


Discrete  Chirp  SAW  ROMs. 

Discrete  chirp  transversal  filters  have  been  used  as  elements  in  CZT 
processor  systems  [4,5].  Such  filters  may  be  regarded  as  acoustic  ROMs  and  used 
to  implement  coherent  frequency  synthesis.  Periodic  impulsing  of  such  a  SAW 
discrete  chirp  filter  will  result  in  periodic  generation  of  the  function 
exp( j 7rn2/N  on  a  carrier,  so  that  two  such  ROMs  operating  at  carrier  frequencies 
f]  and  f2  can  provide  the  necessary  signals  to  generate  a  tone  burst  over  the 
duration  of  the  ROM  outputs.  This  tone  is  obtained  by  delaying  one  signal  with 
respect  to  the  other  and  multiplying  the  two  ROM  outputs: 

p(t-n  At)  cos[2rrfjt  +  rrn^/N]  p(t-n  At)  cos[2  r^t  -  n(n-m)2/N] 

o  0  (2) 

=  p^(t-nAt)  cos[2  7t( f i  +f 2 ) t  +  2rtmn/N  -  nnr/N]  +  (f]-f2)-term, 

where  p(t)  is  the  sampling  window  (e.g.,  p(t)  =  1,  0<t<  At,  and  p(t)  =  0 
otherwise).  The  second  SAW  chirp  ROM  output  is  identical  to  that  of  the  first  one 
except  that  it  is  delayed  by  an  amount  m  At,  and  its  modulation  is  the  complex  conju¬ 
gate  of  that  of  the  first  ROM  (equivalent  to  a  negative  frequency  slope).  Since 

m  can  vary  over  the  interval  0,1 . N-l  in  an  arbitrary  order,  control  over  the 

hopping  sequence  is  flexible. 

Discrete  chirp  ROMs  are  also  useful  since  they  can  be  used  to  generate  a  given 
sinusoid  within  the  band  of  operation  for  an  arbitrary  length  of  time.  This  property 
is  important  for  applications  where  frequency  hopping  is  required  on  a  part-time  bas  ; 
only.  Continuous  SAW  chirp  filters  can  also  be  used  as  ROMs  in  a  similar  sum  or  dif 
ference  frequency  scheme,*  but  care  must  be  taken  to  account  for  possible  discontinuitie: 
where  end-points  of  the  chirp  function  are  encountered. 

Figure  1  illustrates  an  experimental  setup  used  at  NUC  to  perform  frequency 

r 

synthesis  using  SAW  discrete  chirp  ROMs,  where  fi“f2.  The  base-band  output  waveform 
is  shown  in  Figure  2,  and  could  be  filtered  additionally  to  make  the  sampled 
waveform  continuous. 

1 - 

Some  work  in  this  area  has  recently  come  to  our  attention  [11  12], 
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Figure  3  shows  .1  proposed  configuration  for  precise  automatic  control  of  the  hopping 
sequence.  It  makes  use  of  SAW  programmable  diode  multipliers  such  as  have  been 
described  in  [6],  and  incorporates  a  clocking  structure  on  the  same  substrate  such 
as  has  been  described  in  [7,8]. 

Prime  Sequence  ROMs. 

A  third  method  for  ROM  frequency  synthesis  involves  the  use  of  a  serial  access 
memory  in  conjunction  with  a  permuter  so  as  to  achieve  the  equivalent  of  a  random 
access  kOM.  Figure  <1  illustrates  such  a  configuration  for  the  special  case  when  the 
ROM  is  a  SAW  prime-sequence  filter.  Such  a  filter  has  been  designed  and  constructed 
at  NUC  to  implement,  discrete  Fourier  transforms  via  the  prime  transform  algorithm 
[9,10].  this  filter  has  as  its  impulse  response  the  function  exp ( j2 7r p/N) ,  where 
p  =  Rk  mod  N,  k  =  1,2,..., N-l,  and  R  is  a  primitive  root  of  the  prime  number  N. 

Thus,  this  SAW  filter  can  be  used  as  an  acoustic  ROM  to  generate  permuted  values 
of  a  sampled,  complex  sinusoid  on  3  carrier. 

The  attraction  of  this  particular  ordering  of  the  ROM  samples  is  apparent  in  the 
structure  of  the  auxilliary  memory  used  to  operate  the  tap  permuter  (Figure  4).  ihis 
memory  can  be  implemented  <r>  the  form  of  a  cyclically  operated  sequentially- addressed 


binary  ROM  whose  start  ing  position  relative  to  the  pet  iodic  impulsing  of  the  SAW 

acoustic  ROM  determines  the  output  frequency.  N-l  frequencies  can  be  generated,  but 

for  hopping  rates  slower  than  -  1/(N-1)  At  a  complex  sample  value  (1,0)  must  be 

inserted  into  the  synthesized  output  every  N-l  sample'  intervals. 

Such  a  SAW  ROM  synthesis  scheme  might  be  a  reasonable  alternative  to  the 
chirp  method  if  sample  permuters  became  available  whose  performance  exceeded  that 
of  the  chirp  multipliers  in  terms  of  output  accuracy. 


Conclusions. 


The  feasibility  of  using  surface  acoustic  wave  devices  as  principal  elements 
in  ROM  frequency  synthesis  has  been  demonstrated.  In  particular,  the  discrete  chirp 
SAW  device  is  well  suited  because  of  its  ability  to  function  at  the  carrier  of 
interest,  its  inherent  phase  predictability,  its  fast-hop  capability,  and  its  cyclic 

nature  in  the  generation  of  long-duration  single  frequency  waveforms.  Such  SAW 
devices  also  have  the  characteristic  properties  of  relatively  light  weight,  low 
power,  and  small  size,  and  should  be  considered  for  use  as  direct  synthesizers  in 
the  VHP  and  UHF  frequency  regions  and  for  those  applications  where  the  hopping  rate 
required  may  exceed  the  capabilities  of  conventional  digital  techniques. 
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I  XIUMliNTIAL  RESIDUE  CODES 


.1.  M.  A 1  sup  and  J.  M.  Speiser  * 


Abstract  -  This  note  describes  a  class  of  pulse  compression  codes  which 
wore  discovered  by  ex.uuining  the  properties  of  the  prime  transform  algorithm 
[l].  For  each  prime  I’,  phase  hop  sequences  of  length  P-1  can  be  constructed 
whose  periodic  autocorrelation  functions  have  constant  sidelobe  height  of  -1 
relative  to  the  main  peak  value  of  P-1. 


For  each  prime  P,  there  exists  at  least  one  integer  R  such  that  the  resi¬ 
dues  modulo  P  of  R,  R2,  .  .  .  ,  R1*'1  are  all  distinct  and  non-zero,  and  hence 
form  a  permutation  of  1,  2,  ,  P-1  [2,  3].  The  number  of  integers  pos¬ 

sessing  this  property  for  a  given  prime  (P)  is  given  by  the  Euler  0-function, 
<£(P-1)  [2,  4] 1 . 


The  code  corresponding  to  a  selected  integer  (  primitive  root  ')  R  is  the 

»  . 2  o r- 1  , - 

sequence  W  ,  WR  ,  .  .  .  ,  W  ,  where  W  =  exp(j2n/P),  j  =  y-l.  This  sequence 

is  designated  an  exponential  residue  code,  since  individual  terms  are  formed 
by  exponentiation  of  the  residue  of  (j27rRn/P)  modulo  2nj ,  which  exponent  is 
equivalent  to  (j  2zr/P)  (R^iodulo  P) .  This  sequence  can  also  be  recognized  as  a 
particular  scrambled  ordering  of  samples  of  cosine  and  sine  functions  with  sam¬ 
pling  interval  inversely  proportional  to  P. 

The  corresponding  periodic  autocorrelation  function  is 

E1  w(-R")  w(Rk-1)R"=  i  P_1»  k  mod  O'-1)  =  0  (1) 

n=l  n=l  ^-1,  k  otherwise. 

To  prove  the  above  identity,  note  that  when  k  is  between  1  and  P-2  that  Rk  is 
not  equal  to  1,  and  thus  (Rk-1)  is  non-zero.  The  exponents  in  the  sum  therefore 
run  through  the  same  P-1  vlaues  (in  different  order)  as  the  sequence  Rn  itself 
does,  so  that  the  sum  totals  -1. 

To  be  published,  IEEE  Transactions  on  Aerospace  and 
Electronic  Systems.  November  1Q75.  - 
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Since  R,_1  modulo  P  is  congruent  to  1  [2],  the  last  term  in  the  sequence  of 
P-1  elements  is  always  the  element  W,  and  is  always  cyclically  followed  by  the 
element  W*.  Therefore,  the  codes  corresponding  to  different  primitive  roots  are 
distinct,  and  there  are  0  (P-1)  distinct  codes  of  length  P-1. 

Careful  analysis  of  the  Euler  0-function  reveals  that,  for  prime  P>3,  0  (P-1) 
is  always  even.  Furthermore,  each  primitive  root  R  has  an  "associate"  S  belonging 
to  the  set  (1,  2,  .  .  .  ,  P-1)  such  that  the  product  R*S  modulo  P  is  equal  to  one, 
i.e.,  S  *  R-1  [2].  Thus, 

S  =  R'*1  R'1  =  R'-2  '  (2) 


and 


S"=  R~"  =  rp*  ( 1+"/ 


(3) 


That  is,  the  associate  S  is  also  a  primitive  root  which  can  be  used  to  generate  a 
sequence  Ws,  Ws  ,  .  .  .  ,  W*P  1  which  is  equivalent  to  the  sequence  WR  2 ,  W*P  3, 
....  WR,  W.  Thus,  we  can  state  that,  for  P>3,  the  0 (P-1)  primitive  roots  occur 
in  pairs  such  that  one  member  of  a  pair  generates  a  sequence  which  is  the  (shifted) 
reverse  of  that  generated  by  the  other  member. 

An  exponential  residue  code  is  a  shifted  version  of  its  own  complex  conjugate. 
This  can  be  seen  as  follows: 


W 


-(Rrt). 


w 


Rn (P-i) 


wr"i<(',)/,=  ^r" 


i 


(4) 


where  n’  =  n  +  (P-l)/2,  and  the  relationships  lv'P=l;  (P-1)2  modulo  P  -  1 ;  and  Rp~‘ 
modulo  P  =  1  have  been  utilized.  The  corresponding  correlation  relation  also  holds: 

P-1  nj"+S  lP-1,  k  mod  (P-1)  =  (P-l)/2 

Z  W(R  )  W<R  >  =  <  (5) 

n=l  /  -1,  k  otherwise. 

The  autocorrelation  function  of  the  exponential  residue  phase  hop  codes  may 
also  be  derived  by  considering  Rader's  prime  transform  algorithm  [1]  when  the  sig¬ 
nal  to  be  analyzed  is  one  of  the  basis  vectors  of  the  discrete  Fourier  transform, 
and  hence  the  output  of  the  prime  transform  is  a  Kronecker  delta.  Other  pulse 
compression  codes  are  known  which  use  the  analog  to  a  logarithm  in  modulo  P 
arithmetic  [5,  6].  Frequency-hop  codes  of  a  type  related  to  exponential  residue 
phase  hop  codes  are  also  of  interest  and  have  been  examined  by  Cooper  and  Yates  [7J. 


‘From  [2],  when  a  number  A  is  written  in  terms  of  its  prime  factors  a,  b,  c,  etc: 
then  4(4)  is  given  by 


A  =  aa  b^  cY  etc, 


etc- 


(6) 

(7) 
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